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Abstract 


Direct numerical simulations have been used to examine the effect of the initial distur- 
bance field on the development of three-dimensionality and the transition to turbu- 
lence in the incompressible plane wake. The simulations were performed using a new 
numerical method for solving the time-dependent, three-dimensional, incompressible 
Navier-Stokes equations in flows with one infinite and two periodic directions. The 
method uses standard Fast Fourier Transforms and is applicable to cases where the 
vorticity field is compact in the infinite direction. Initial disturbances fields exam- 
ined were combinations of two-dimensional waves and symmetric pairs of 60° oblique 
waves at the fundamental, subharmonic, and sub-subharmonic wavelengths. 

The results of these simulations indicate that the presence of 60° disturbances 
at the subharmonic streamwise wavelength results in the development of strong co- 
herent three-dimensional structures. The resulting strong three-dimensional rate- 
of-strain triggers the growth of intense fine scale motions. Wakes initiated with 
60° disturbances at the fundamental streamwise wavelength develop weak coherent 
streamwise structures, and do not develop significant fine scale motions, even at high 
Reynolds numbers. The wakes which develop strong three-dimensional structures ex- 
hibit growth rates on par with experimentally observed turbulent plane wakes. Wakes 
which develop only weak three-dimensional structures exhibit significantly lower late 
time growth rates. 

Preliminary studies of wakes initiated with an oblique fundamental and a two- 
dimensional subharmonic, which develop asymmetric coherent oblique structures at 
the subharmonic wavelength, indicate that significant fine scale motions only develop 
if the resulting oblique structures are above an angle of approximately 45°. 
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Chapter 1 


Introduction 


1.1 Background 

This section presents a review of results from previous theoretical, experimental, and 
computational studies of incompressible plane wakes. A brief review of numerical 
simulations of free shear flow’s is also presented. 

1.1.1 The Incompressible Wake 

Significant progress has been made in understanding the primary stages of transition 
in incompressible wakes behind slender bodies. There exists substantial theoretical, 
experimental, and computational work describing the initial development of the initial 
wake instability: the growth of two-dimensional Kelvin-Helmholz w’aves into the well 
known Karman vortex street. 

Sato & Kuriki [39] performed a now classic set of experiments on both natural and 
forced wakes behind a thin flat plate. They identified three distinct developmental 
regions in the wake: the linear region (also known as the Kelvin-Helmholz instability 
region), where small amplitude disturbances grow exponentially; the nonlinear growth 
region, w’here the fundamental two-dimensional mode saturates and the wake rolls up 
into the Karman vortex street; and the three-dimensional region, w’here strong three- 
dimensional motions appear. They found that growth rates and mode shapes of the 
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disturbances in the linear region were well predicted by linear stability analysis of 
a Gaussian mean base flow. Past the short linear region, however, they observed 
that the evolution of the wake deviates substantially from the predictions of linear 
theory. The amplitude of the two-dimensional fundamental disturbance saturated, 
and then decreased, higher harmonics of the two-dimensional fundamental appeared 
in the w r ake, and the mean velocity and wake width changed at a rate which was far 
more rapid than could be accounted for by linear theory. 

Ko, Kubota, k Lees [22] performed a two-dimensional finite amplitude, single 
frequency disturbance analysis of the plane wake to explain the results reported bj 
Sato k Kuriki for the nonlinear region. By applying an integral method to a bound- 
ary laver approximation of the Navier-Stokes equations, they were able to study the 
energv balance between the mean flow and the finite amplitude disturbances. The\ 
related the observed behavior of the wake disturbance amplitude and the rapid varia- 
tions in the mean flow to energy transfer between the disturbance and the mean flow 
via the Reynolds stresses. Their analysis also emphasized the importance of binary 
interactions between disturbance modes, which are the source of higher harmonics in 
the wake. 

These studies, plus others such as those by Sato k Saito [40], who extended the 
work of Sato k Kuriki [39] by examining the effect of multiple frequency forcing, and 
Mattingly k Criminale [26], who developed a disturbance theory which included the 
effects of the non-parallel nature of the near-field wake, give a fairly comprehensive 
picture of the early stages of wake instability. Understanding of two-dimensional near 
wakes behind bluff bodies much more difficult since the initial instability is dominated 
by a complex mixture of large amplitude motions including flow separation dynamics, 
rollup of the shed shear layers, and reversed flow in the near wake (c.f. the studies of 
Karniadakis & Triantafyllou [20, 21]). It is reasonable to assume, however, that far 
downstream, away from the vicinity of the body, the dynamics of the fully developed 
blunt body wake should be similar to those for the wakes of slender bodies. Only the 
initial conditions should differ. 

In contrast to the fairly complete understanding of the initial stability characteris- 
tics and early development of the plane wake, there still remain significant unresolved 
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issues related to the stability characteristics of the two-dimensional non-linear wake. 
It is these secondary stages of instability that lead to the appearance of strong three- 
dimensional motions in the wake. 

Early experiments on the three-dimensional structure of the far wake were per- 
formed by Townsend [45]. These experiments were later extended by Grant [17], who 
took long-time-averaged velocity correlations in the far wake of a cylinder. Coherent 
three-dimensional structures, later referred to as “double-roller eddies” by Townsend 
[46], were inferred from these measurements. These structures were described as pairs 
of curved, counter-rotating vortices oriented perpendicular to the plane of the wake. 
Roshko [37] suggested that the time-averaged data was actually due to a superposi- 
tion of vortex loops, formed by the distortion of spanwise vortices from opposite sides 
of the wake. Mumford [29] later used a pattern recognition technique to study the 
double-roller eddies, and concluded that they were often confined to one side of the 
wake or the other, and that they tended to appear in groups. 

A set of experiments by Cimbala, Nagib, k Roshko [13] were among the first 
in which the existence of spanw'ise periodic streamwise structures in the planar far 
wake of a bluff body was documented. They identified the structures as hairpin 
vortices produced by a parametric resonance between the two-dimensional and oblique 
subharmonic disturbances. They regarded this secondary instability as similar to one 
studied by Pierrehumbert k Widnall [32] in plane mixing layers. 

Flemming [16] performed an analysis of the secondary instability of the plane 
wake. Taking as a base flow a Gaussian mean velocity field with a two-dimension- 
al fundamental Orr-Sommerfeld mode superimposed, Flemming obtained a Hill-type 
system of equations for the stability of the far wake. Numerical studies of the stability 
equations indicated that, for sufficiently high fundamental (Orr-Sommerfeld) mode 
amplitudes, pairs of oblique waves at the subharmonic streamwise wavelength were 
unstable. The angle of the most unstable disturbance was found to depend on the 
wake Reynolds number, but unstable disturbances were found to exist at angles of 
between 45° and 70° with respect to the spanwise direction. 

Corke, Krull, k Ghassemi [14], using the results of Flemming, performed a study 
of the mechanisms for the secondary growth of three dimensional modes in the far 
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wake of an airfoil. They focused on a parametric resonance mechanism between the 
fundamental two dimensional Karman instability mode and pairs of phase locked 
oblique disturbances at the subharmonic wavelength that were oriented 60° and 73° 
with respect to the spanwise direction. These resonances were expected from the form 
of the stability equations developed by Flemming. They found that, at least under 
certain conditions, a resonance develops where the two dimensional fundamental mode 
and the subharmonic oblique modes exchange energy over several long period cycles. 

Williamson [52] and Williamson k Prasad [53, 54, 55] suggest an alternate mech- 
anism for the development of strong three-dimensional motions. Recent experiments 
they performed indicate that the oblique waves observed in the far wake of a cylinder 
by Cimbala et a 1. are due to an interaction between oblique shedding waves produced 
by the wake generator and two-dimensional subharmonic waves which arise from the 
inherent hydrodynamic instability of the mean flow in the far wake. 

These differing theories on the source of three-dimensional motions — one which 
ascribes the three-dimensional motions in the far wake to the growth of pairs of highly 
oblique subharmonic waves which exist from early in the flow, and one which ascribes 
the three-dimensional motions in the far wake as being composed of oblique waves 
generated from an interaction between asymmetric fundamental vortex shedding and 
long wavelength two-dimensional motions — reinforce the need to examine the effect 
of initial conditions on the development of the wake. Both theories suggest the 
importance of studying the dynamics of highly oblique disturbances in plane wakes. 

Wygnanski, Champagne, k Marasli [58] conducted an experimental study of small 
deficit turbulent wakes created using a variety of wake generators. The generators 
were carefully chosen to have the same drag, and therefore to create wakes with 
the same momentum thickness and momentum thickness Reynolds number. They 
found that the normalized characteristic velocity and length scales as well as the 
normalized longitudinal turbulence intensity depended on the generator used, and 
hence depended on the initial conditions created by that generator. The shape of 
the mean velocity profile, however, was found to be independent of initial conditions. 
They attributed the lack of universality in part to the interaction between the sinuous 
(cross-stream component of velocity antisymmetric about the wake centerline) and 
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varicose (cross-stream component of velocity symmetric about the wake centerline) 
instability modes in the far wake. A common interpretation of classical similarity 
theory, w’hich characterizes the turbulent wake through a single parameter the 
momentum deficit (see appendix A) — holds that the turbulent wake should have 
a universal growth rate, independent of initial conditions. In fact, classical theory 
addresses only the late time growth laws for the turbulent wake, not the actual rates. 
The implication of the results of Wygnanski et a I. is that the rate of development of 
the turbulent far field of the plane wake is indeed dependent on parameters related 
to the initial conditions, and therefore not universal. 

In addition to the primarily experimental and theoretical studies discussed above, 
there have been a large number of computational studies of free shear flows. 

Mixing layer simulations by Riley, Mourad, Moser k Rogers [33] indicates that 
the large scale structures that develop in an incompressible mixing layer are strongly 
dependent on the phase angle between a two-dimensional fundamental disturbance 
and a pair of oblique disturbances. They found that the existence of strong vorticity in 
the region between spanwise rollers at late times required the presence of streamwise 
vorticity in the same region from early on in the development of the layer. They 
determined that the intense streamwise vorticity was produced by stretching of the 
early vorticity by the strain field induced by the large structures. 

Moser k Rogers [27] and Rogers k Moser [35, 36] conducted a comprehensive 
numerical study of a temporally evolving incompressible plane mixing layer started 
from “clean’’ initial conditions (the initial conditions consisted of a mean flow plus 
a small number of low wavenumber disturbances). They found that most of the sets 
of initial conditions they studied led to the development of very intense streamwise 
vortical structures. Those mixing layers that did not develop these strong streamwise 
structures during the initial period of growth took much longer to develop any signif- 
icant three-dimensionality. Because of this they concluded that the development of 
the strong streamwise structures were a key step in the development of three-dimen- 
sionality and the eventual transition to turbulence. 
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Chen, Cantwell, k Mansour [9, 10] carried out a direct numerical simulation of 
a temporally evolving compressible plane wake. They found that linear theory ac- 
curately predicted the early growth of the plane wake for freestream Mach numbers 
between M =0.01 and M = 3.0. They also found that the development of three-dim- 
ensionality in the compressible wake was significantly affected by the relative phase 
of the initial disturbance. 

Maekawa, Mansour, k Buell [24] performed direct numerical simulations of a set 
of two-dimensional spatially evolving incompressible plane wakes. They found that 
wakes initiated with a two-dimensional fundamental and a two-dimensional subhar- 
monic disturbance initially form a Karman vortex street at the fundamental wave- 
length. Once the fundamental has saturated, the subharmonic disturbance begins 
to become significant and the vortices in the vortex street combine, forming pairs 
of both like and opposite signed vortices. A wake forced with a combination of a 
fundamental and random noise showed similar behavior. 

Moser &: Rogers [28] performed a direct numerical simulation of a pair of tempo- 
rally evolving incompressible plane wakes started from pairs of temporally evolving 
turbulent boundary layers which had been previously computed by Spalart [43]. The 
first wake was initiated with only the computed boundary layers. The second was 
initiated with modified boundary layers which had all the two-dimensional modes in- 
creased by a factor of 20 in an attempt to simulate the receptivity of the plane wake 
to two-dimensional disturbances which normally would occur at the wake generator. 
The high amplification of the two-dimensional modes, which amounted to a 13 fold in- 
crease in total disturbance energy, was found to be necessary to spur the development 
of the large scale two-dimensional structures that were expected to appear. 

The first turbulent wake case computed by Moser k Rogers showed eventual self- 
similar (t 1 / 2 ) growth and energy spectra with a short region of k ~ 5/3 slope. However, 
the late time growth rate was found to be well below the range of rates measured 
by Wygnanski et a 1. [58]. The second wake, with the enhanced two-dimensional 

disturbance field, also showed a region of k ~ 5 ^ 3 spectra, but never developed the 
expected self-similar t growth pattern for any extended period of time. The forced 
wake also had a growth rate well above the Wygnanski range. 
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1.1.2 Topology of Fine Scale Motions 

One of the major unsolved problems of fluid mechanics is how to model turbulent 
flows. Model development has been hindered by the fact that different types of flows 
require different models since the large scale features of the turbulence are highly flow 
dependent. An attempt to get around this problem is the technique of large eddy 
simulation, which splits the flow into the mean flow and large scale turbulence, which 
is simulated numerically, and the fine scale turbulence, which is modeled. The basic 
premise of the technique is that the fine scale turbulence has features which are flow 
independent and therefore more amenable to modeling than the larger scales. 

The typical approach to developing a model for the fine scale turbulence is to 
assume statistical isotropy of the turbulent motions at high wavenumbers there- 
fore relying on the assumption that the spectral characteristics of the turbulence are 
universal. A number of recent studies of the topology of the fine scale velocity fields 
of a variety of turbulent flows have revealed another potential path to developing a 
turbulence model for use in large eddy simulations. They have found what appear to 
be universal features in the geometric properties of fine scale turbulence in physical 
space. The existence of such universal features could potentially lead to models for 
fine scale turbulence based on the physical (local) rather than the spectral (global) 
properties of the turbulence. 

Ashurst, Kerstein, Kerr k Gibson [2] studied direct numerical simulations of in- 
compressible forced isotropic turbulence and homogeneous sheared turbulence. They 
found that the intermediate principal strain-rate tended to be positive throughout 
the flow. Furthermore, data conditioned on high levels of local dissipation of kinetic 
energy had a uniformly positive intermediate strain-rate, with strain rates in the ra- 
tio of approximately 3:1: —4. They also found that the highly dissipating motions 
tended to have the vorticity vector aligned with the intermediate principal strain-rate 
direction. 

More recent studies of incompressible forced isotropic turbulence by Vincent k 
Meneguzzi [50, 51] and Ruetsch k Maxey [38] indicate that the small scale structures 
take the form of vortex tubes. The highest rates-of-strain, and therefore highest rates 
of dissipation, were found to occur in the vicinity of, but outside the cores of. these 
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vortex tubes. Vincent k Meneguzzi again found that the vorticity in these high dis- 
sipation regions was aligned with the intermediate principal strain-rate. Furthermore 
they found that this result held before the vortex tubes had developed. 

Chen et a I. [11] studied a set of direct numerical simulations of transitional com- 
pressible and incompressible mixing layers. They used a new data display method 
based on the classification of local flow topologies using the techniques outlined by 
Chong, Perry k Cantwell [12] (see appendix D). The method allows for simple and 
straightforward study of global trends in the topology of the fine scales and correla- 
tions between physical and topological features. As in the other studies, they found 
that the intermediate principal rate-of-strain tended to be positive, and that the trend 
became stronger as the data was conditioned on higher rate of dissipation. They also 
found that regions of high dissipation tended to be associated with similarly high 
enstrophy density. Sondergaard et a 1. [41], Soria, Sondergaard k Cantwell [42] and 
Blackburn, Mansour k Cantwell [3] extended the study of Chen et a 1. to include 
data from simulations of compressible and incompressible wakes, a turbulent incom- 
pressible mixing layer, and incompressible channel flow respectively. Again the same 
general topological features were observed. 

Finally, Tsinober, Kit k Dracos [47] performed an experimental study of the align- 
ment of strain and vorticity in both grid-generated and boundary layer turbulence. 
They observed a tendency for the vorticity vector to align with the intermediate prin- 
cipal rate-of-strain in agreement with the previous studies of numerical simulations. 

Attempts have been made to explain these observations. Jimenez [19] suggested 
a kinematic model for the alignment of the vorticity and strain using the stretched 
Burgers’ vortex as an example. Though the model described a vortical flow in which 
the observed alignment occurred, there tvas no attempt to explain the evolution of 
such structures in a real flow. Cantwell [5] studied a restricted Euler equation, first 
studied by Vieillefosse [48, 49], in which viscous terms and mixed second derivatives of 
the pressure had been dropped. The resulting closed-form solution for the evolution 
of the velocity gradients reproduced the tendency for the strain-rates to evolve to a 
state with a positive intermediate principal rate-of-strain. Cantwell [6] later developed 
an intermediate asymptotic model for the case where the viscous terms and mixed 
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derivatives of the pressure were non-zero. This model helped to explain the structure 
of the invariant pdfs in the mixing layer data examined by Soria et a I. 


1.1.3 Goals of This Study 

In light of the studies discussed above, the following questions arise: 

• Is the parametric resonance model proposed by Flemming and Corke et a 1. for 
the development of three-dimensionality in the far wake an appropriate one? If 
so: 

• How do the initial conditions affect the development of three-dimensionality in 
the incompressible plane wake? 

• How do the initial conditions affect the mean flow and structure of the turbu- 
lence in the far wake? 

• How do the initial conditions affect the growth rate and mean properties of the 
turbulence in the far wake? 

The intent of this study is to begin to address these questions. 

1.1.4 Direct Numerical Simulation of Free Shear Flows 

The basic tool used in this study is direct numerical simulation. Here the term “di- 
rect” indicates that there is no attempt to model unresolved scales in the simulation. 
All of the scales in the flow which contain significant amounts of energy are numeri- 
cally resolved and evolve as solutions to the full Navier-Stokes equations. 

The major shortcoming of numerical simulation is the limitation on resolution. A 
numerical simulation is restricted by the size and speed of the computational hardware 
used. For a given simulation on a given machine, there is a fixed range of scales which 
can be reasonably computed. For the time being at least, this means that flows studied 
bv numerical methods in general, and direct numerical methods in particular, are 
either limited to relatively small flow domains or to Reynolds numbers that are quite 
low compared to those obtainable in laboratory experiments. 
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The size of the flow domain which can be represented numerically is also limited 
by the computational resources available. This limitation can have an impact on the 
development of the computed flow. While this constraint also appears in laboratory 
experiments (w'here walls, boundary layers, and limited facility length can affect the 
flow), it is typically less severe in experiments. Because of these facts, care must be 
taken when attempting to generalize the results of numerical studies. 

To maximize the range of computed flow scales, a temporal formulation has been 
used in the simulations performed for this study. A temporal simulation may be 
thought of as approximating the evolution of a representative set of structures in 
the physically realizable spatially evolving flow as they convect downstream. In the 
case of the wake, a temporal formulation approximates the view an observer that was 
convecting downstream with the freestream (or alternately a fixed observer that has 
been passed bv the wake generator) would have of the evolution of the flow structures. 
In a temporal formulation, the roles of time and the downstream coordinate direction 
are swapped with respect to the corresponding spatial formulation. Time in the 
temporal formulations becomes the measure of the level of development of the flow, 
in place of downstream distance for the spatial formulation. Variations in the flow 
at different streamwise coordinates at a fixed time in the temporal formulation are 
analogous to variations with time at a fixed point in the spatial formulation. 

Using a temporal rather than a spatial formulation has the advantage of allowing 
the resolution of smaller scales for a given set of computational resources. A spatial 
formulation requires the resolution of both a large flow domain and the fine scale 
motions. This represents a potentially very wide range of scales, with correspondingly 
large computational requirements. Using a temporal formulation allows the range of 
resolved scales to be reduced by restricting the largest resolved scales in the flow to 
at most a few representative large structures instead of the full domain. By reducing 
the computational requirements associated w'ith the large scales, more of the available 
resources may be focused on resolving the small scales. 

One assumption inherent in any temporal formulation is that the streamwise rate 
of change of a spatially evolving flow is small at the scale of the structures being 
studied. When this assumption holds, the mean flow may be approximated as being 
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locally parallel without significantly affecting the development of the structures of 
interest. This assumption is generally a very poor one very near the origin in any 
spatially developing free shear flow. It can, however, be quite reasonable away from 
that region if the growth of the far field flow is sufficiently slow. For wakes, in 
which the far downstream flow grows as (where x is the downstream distance) 
the assumption of parallel mean flow is quite good. For mixing layers, in which the 
far downstream flow grows like x, the assumption of parallel mean flow is less valid. 
Under the proper conditions, however, the change in width of the mixing layer over 
the length of the structures being studied can be relatively small, making the temporal 
approximation useful. 

Temporal formulations also inherently ignore the details of how disturbances are 
initially created in the flow. Since here we are interested in how specific disturbance 
modes (disturbances with specific wavelengths, angles, and phases) affect the devel- 
opment of the far wake, the exact details of how those specific modes are generated in 
a spatially evolving wake are not addressed. All of the geometric details of the wake 
generator are subsumed in the choice of disturbance modes and disturbance phases. 

Finally, temporal formulations (and to a lesser extent spatial formulations), by 
limiting the computational domain to a few large structures in the flow, can confine 
the development of the flow in the streamwise and cross-stream directions, possibly 
quite severely. This artificial confinement can inhibit or prohibit the development 
large scale motions which might normally exist, and become dynamically significant, 
in a physical flow. The effects of the finite computational domain must be taken into 
account when interpreting the computational results, particularly in the late stages 
of development of the flow.. 

Numerical simulations in general have some special advantages in answering the 
type of questions asked at the end of section 1.1.1 above. The initial conditions can be 
very precisely controlled and the results present a complete description of the entire 
flow field, including physical variables that would normally be very hard to measure. 
This allows for unambiguous connections to be made between initial conditions and 
developed structures. In addition, direct numerical simulations give access to all the 
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physical quantities in the flow at a given instant. This allows study of quantities 
which are not normally available from a laboratory experiment. 

The numerical technique used here to perform the simulations is known as a 
“pseudo-spectral” method. In a “spectral” method, the dependent variables are ex- 
panded as a sum of (usually orthogonal) basis functions. This allows the governing 
partial differential equations for the physical variables in the problem to be converted 
into a set of ordinary differential equations for the time evolution of the coefficients 
of the basis functions in the approximating sums. The primary benefit of such an 
approach is that, by expressing spatial derivatives of the basis functions in terms 
of the basis functions, calculation of spatial derivatives is greatly simplified. It also 
allows derivatives to be calculated with “spectral’ accuracy. This means that the 
error in the representation of derivatives goes to zero exponentially as the number of 
functions in the basis set goes to infinity. Other approaches, such as finite difference 
schemes, typically have errors that go to zero algebraically. 

Another benefit of spectral methods is that they lead to algorithms that are simple 
to implement on parallel processing computers (machines that are designed to use 
multiple interconnected microprocessors to work different parts of the same problem 
simultaneously). Larger and more complex simulations can be carried out on these 
parallel machines. 

In what is known as a “fully” spectral method, any nonlinear terms in the gov- 
erning equations are computed using convolution integrals involving the coefficients 
of the expansions. This is a computationally intensive process. A “pseudo”-spectral 
method makes use of fast transforms to convert between the physical representation of 
the dependent variables and the basis function expansions. This allows any nonlinear 
terms to be calculated in physical space using simple multiplication then re-expanded 
in terms of the basis functions. For large problems, the use of fast Fourier trans- 
forms makes pseudo-spectral methods significantly more efficient than fully spectral 
methods. 
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1.2 Outline of Present Work 

The intent of the present work is to examine the effect of the choice of initial conditions 
on the development of the incompressible plane wake. Of particular interest is how 
the initial conditions relate to the development of three-dimensionality and eventual 
transition to turbulence in the wake. Also of interest is whether the structure of the 
turbulence which develops is independent of initial conditions as predicted by the 
usual one parameter similarity analysis of the far wake. 

All of the flows that will be described here were started from a laminar base 
flow upon which was superimposed a small number of disturbance modes at very 
specific wavelengths. The disturbance modes were the most unstable eigenfunctions 
as predicted by linear stability theory. Understanding these very simple, clean , 
wake flows should allow for a better understanding of wakes started with more realistic 
initial conditions, which will contain uncontrolled disturbances. 

Chapter 2 describes the numerical methodology used to perform the direct numer- 
ical simulations used in this study. A new pseudo-spectral algorithm for simulating 
planar shear flows with periodic freestream boundary conditions is described and 
tested. The method uses Fourier transforms in all three spatial directions to solve for 
the flow in a finite, time- varying, computational domain using velocities which are 
matched at the domain boundaries to known analytic solutions. 

Chapter 3 presents results from a series of two-dimensional direct numerical simu- 
lations. The effect of the choice and phasing of two-dimensional disturbance modes is 
described. The effect of flow Reynolds number on the development is also examined. 

Chapter 4 presents results from a wide series of three-dimensional direct numer- 
ical simulations of the plane incompressible wake. The effects of disturbance mode, 
relative disturbance phasing, and flow Reynolds number on the wake development 
are discussed. 

Chapter 5 examines the evolution of the topology of the small scale motions for a 
selected subset of the three-dimensional simulations discussed in chapter 4. 

Chapter 6 presents the major conclusions of this work and outlines some recom- 
mendations for future work. 
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Appendix A presents a brief outline of classical similarity theory as applied to 
incompressible plane wakes, as well a definitions of some of the turbulence measures 
used. 

Appendix B presents a review of basic linear stability theory and the method used 
to generate the disturbance eigenfunctions for the computed flows. 

Appendix C presents an overview of numerical aliasing and a description of the 
dealiasing algorithm used. 

Appendix D describes in greater detail the topological methods used in chapter 5 
to study the structure of the fine scale motions in the computed flows. 

Appendix E presents a full set of invariant space pdfs for selected wake simulations. 

Appendix F presents mean turbulence statistics for selected wake simulations. 

Appendix G presents preliminary results from a set of wake simulations initiated 
with combinations of oblique fundamentals at various angles and two-dimensional 
subharmonic disturbances. This approximates the initial conditions in the experi- 
ments of Williamson and Williamson & Prasad. 

Appendix H gives a listing of all the computations performed for this study along 
with the values of the relevant flow parameters for each simulation. 


1.3 Summary of Results 

1.3.1 The Two-dimensional Plane Wake 

Two-dimensional simulations of the temporally evolving plane wake initiated with 
combinations of disturbance eigenfunctions at the fundamental, subharmonic, and 
sub-subharmonic wavelength indicate the following: 

• The presence of a subharmonic disturbance causes the initial Karman rollers 
to amalgamate and/or pair, depending on the phase relative to the fundamen- 
tal. The presence of a sub-subharmonic disturbance has minimal effect on the 
structure or growth of the wake. 

• Increasing the Reynolds number increases the intensity of the large scale struc- 
tures in the flow at a given development time without having a major effect 
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on their overall shape. Increasing Reynolds number significantly increases the 
number, intensity, and duration of small scale features. 

• The passive scalar and enstrophy density fields track each other well. Major 
differences between the enstrophy density and scalar fields only appear in regions 
of the flow where there has been significant cancellation of vorticity of opposite 
sign. 


1.3.2 The Three-dimensional Plane Wake 

Physical Space 

Temporal simulations of the three-dimensional plane wake at various Reynolds num- 
bers, initiated with combinations of two-dimensional and three-dimensional distur- 
bance eigenfunctions at the fundamental, subharmonic, and sub-subharmonic wave- 
length at various phases indicate the following: 

• The mechanism proposed by Flemming and Corke et a I. is a legitimate route 
for the development of three-dimensional motions in the far wake. 

• Wakes with two-dimensional and oblique disturbances at only the fundamental 
wavelength do not produce any significant three-dimensionality in the far wake. 
The addition of a two-dimensional subharmonic disturbance produces coherent 
three-dimensional structures of only moderate strength. 

• The presence of an oblique disturbance at the subharmonic wavelength results 
in the development of very strong three-dimensional structures, independent of 
the presence or absence of disturbances other than the two-dimensional funda- 
mental. 

• The primary effect of the phase of the oblique disturbances is to determine 
which side of the wake develops the dominant three-dimensional structures. The 
phase of the two-dimensional subharmonic disturbances has a significant effect 
on the development of streamwise structures in the wake only in the absence 
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of subharmonic oblique disturbances (when the w r ake dynamics are primarily 
two-dimensional). 

• As the wake Reynolds number is increased, the strength of the three-dimen- 
sional structures increases. At the highest Reynolds numbers simulated, where 
there is a significant volume of intense vorticity spread through the wake, the 
highest enstrophy density regions appear as coherent three-dimensional struc- 
tures. These structures take the form of elongated vortex tubes with lengths on 
the order of the wake width and diameters on the order of the turbulent scales. 

• Preliminary studies indicate that the mechanism proposed by Williamson k. 
Prasad for the development of three-dimensional motions in the far wake is 
also legitimate so long as the shedding angle of the oblique fundamental is 
sufficiently high. Weakly oblique shedding does not appear to result in the 
strong streamwise structures that are necessary for the development of fine 
scales. 

Invariant Space 

The simulation results have been used to study the invariants of the velocity gra- 
dient tensor. Topological analysis of the fine scale, high gradient, motions in the 
incompressible wakes revealed the following: 

• The wakes with a three-dimensional subharmonic have both a greater quantity 
of and more intense high gradient motions. 

• The characteristic shapes of the joint probability density functions for invariants 
of the incompressible plane wake are similar to those observed in the other three- 
dimensional flows discussed in section 1.1.2: 

— Joint pdfs of the second and third invariants of the velocity gradient tensor 
have a characteristic “skewed teardrop” shape, with high gradient motions 
tending to be of topological types stable vortex /stretching and unstable 
node/ saddle/ saddle. 
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- Joint pdfs of the second and third invariants of the rate-of-strain tensor 
indicate that the most dissipative motions are associated exclusively with 
an unstable node/saddle/saddle type strain topology. More moderately 
dissipating motions, which account for the majority of the integrated dis- 
sipation in the flow, are also very strongly associated that strain topology. 

- Joint pdfs of the enstrophy density and vortex stretching indicate that 
highly rotational motions occur in regions where the vortex stretching is 
positive, even at late times when the intensity of all gradients are de- 
creasing in the wake. In addition, regions with moderate to high rates of 
dissipation tend to have the vorticity vector aligned with the intermediate 
principal rate-of-strain direction. 

• Increasing the wake Reynolds number increased the intensity of the gradients 
while preserving the shape of the pdfs in invariant space. In effect, the shape of 
the pdfs are Reynolds number invariant. In addition, changes in the phases of 
the initial disturbances have minimal effect on the shape of the invariant space 
pdfs. 



Chapter 2 

Numerical Methodology 


2.1 Introduction 

The challenge of calculating free shear flows using spectral methods has been ap- 
proached using a variety of numerical schemes. The problems of interest are turbulent 
flow's which are periodic in two directions and have vorticity which is “compact” (of 
finite extent) in the third direction. An early computation of such a case was per- 
formed by Orszag k Pao [30] who simulated a temporally developing momentumless 
wake using a pseudo-spectral method. They approximated the infinite direction in 
the flow by truncating to a finite domain and using sine and cosine transforms in the 
inhomogeneous direction, effectively imposing free slip conditions at the non-periodic 
boundaries of the finite computational box. A similar approach using a three-dim- 
ensional vorticity stream function formulation was taken by Mansour, Ferziger, k 
Reynolds [25] to compute a time-developing turbulent mixing layer using a large 
eddy simulation technique. The disadvantage of approaches such as these is that 
they fail to accurately treat the irrotational field in the infinite direction by forcing 
the flow to be parallel some finite distance into the freestream regions. 

Cain, Reynolds, k Ferziger [4] analyzed the method of Orszag k Pao [30] and 
found that the size of the domain in the inhomogeneous direction can influence the 
computational results. To circumvent the restriction of a finite size domain, they 
introduced a cotangent mapping in the infinite direction which allowed the use of 
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Fourier spectral methods in the doubly infinite domain. This effectively moved the no- 
stress boundaries out to a very large distance from the rotational region, minimizing 
the error inherent in such an artificial boundary condition. This sort of approach 
has become the mainstay of numerical simulations of flows with infinite domains (i.e. 
Chen, Cantwell & Mansour [10], Laurien &; Kleiser [23]). The disadvantage of this 
scheme is that it sacrifices some of the simplicities of using the Fourier transform (such 
as differentiation in physical space being represented by simple multiplication by a 
wave vector in Fourier space), resulting in greater coding complexity. Computational 
resolution is also wasted on regions of the flow which are free of vorticity. 

Spalart, Moser, & Rogers [44] approached the problem by using a set of basis 
functions for their expansions that are defined on the semi-infinite interval. They used 
a set of Jacobi polynomials in a mapped variable to represent the vortical region, and 
slowly decaying exponential “extra” functions to accurately represent the potential 
flow far from the vortical region. This combination of quickly decaying and slowly 
decaying basis functions allowed them to achieve good accuracy and good (though 
not spectral) convergence. The drawback is that the use of Jacobi polynomials is 
numerically expensive as each Jacobi transform must be accomplished by quadrature. 
In addition, this approach, as well as the Cain mapping approach, tends to concentrate 
resolution very near the centerline, at the expense of resolution a small but finite 
distance away. This approach is well suited to mixing layers, but is inconvenient for 
flows for which high gradients and small scales develop away from the centerline as 
occurs in the wake. 

The goal of the method described here is to avoid the drawbacks of the previous 
methods. It is based on an algorithm similar to one presented by Corral & Jimenez 
[15]. Fourier transform techniques, for which fast numerical transforms exist, are used 
to solve problems which have one infinite and two periodic directions without hav- 
ing to resort to nonlinear mapping of the infinite direction to a finite domain. This 
preserves all of the benefits of using Fourier transforms (differentiation accomplished 
by multiplication by a wave number, integration accomplished by division by a wave 
number, interpolation accomplished by multiplication by a phase factor, resolution 
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Figure 2.1: Numerical domains 

changes accomplished by adding or truncating zeros in wave space) while still accu- 
rately representing the boundary conditions. It also allows for uniform resolution of 
the domain of interest. 

2.2 Approach 

The basic approach is to divide the flow into three domains in the non-periodic 
direction as shown in Figure 2.1. Domain I extends from — oo to Vi(f), domain II 
from Y A (t) to V^f), and domain III from V^f) to +oo. The boundaries Vi(f) and 
are chosen such that domain II contains all the vorticity in the flow, and domains I 
and III are vorticity free. This choice of yi(f) and 1^(0 requires that the vorticity 
magnitude and vorticity gradient at the top and bottom of domain II be zero (or 
in practice very small). This is the case for many flows of interest, particularly at 
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moderate to high Reynolds number where the interface between rotational fluid and 
irrotational fluid is sharp. Choosing Vj(i) and in this way allows the vorticity in 
domain II to be treated as fully periodic. The vorticity equations governing the flow 
may then be solved with a pseudo-spectral technique which uses standard complex 
Fourier transforms in all three directions. 

The inherent non-periodic nature of the flow in the cross stream direction only 
enters into the equations of motion through the nonlinear term of the momentum 
equation, which involve the velocity. The vorticity, u > : , is not directly affected by the 
images of the flow created by artificially introducing periodicity in the non-periodic 
direction. The vorticity is effectively zero at the nonperiodic edges of the box and 
can accurately be expanded using periodic functions. The velocity, Uj, which is a 
solution of a Poisson equation involving the vorticity as a source term, is affected 
by the vorticity images and must be corrected to remove the effect of the artificial 
periodicity. This is accomplished by adding an incompressible, irrotational component 
to the velocity field in domain II which matches it to analytic asymptotic solutions 
for the velocity in domains I and III. 

Note: In the following discussion, all quantities have been normalized by the 
initial flow halfwidth b 0 and the freestream velocity U 0 as follows: 


_ _»£ £ = J_M # 

Jj ~bo U] ~U 0 b 0 p Ug{p) 

"° (0 n 

V b 0 U 0 K b 0 U 0 ~ Pr ( ) 

where ()° are the unnormalized quantities. The initial halfwidth, bo is defined as half 

the width of the initial mean velocity profile at half the maximum mean defect velocity 

(see figure 2.2). The mean profile is generated by averaging over X\-Xz planes. 
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Figure 2.2: Mean flow parameters 


2.3 Governing Equations 

2.3.1 Vorticity Form of the Navier-Stokes Equations 

The incompressible, uniform density, Navier-Stokes equations are 

ujj = 0 (2.2) 

Uj, t + u k u hk + — = vu hkk (2.3) 

P 

where u 3 is the velocity in the x 3 direction, p is the pressure, and p is the (con- 
stant) density. Here, and throughout the remainder of this dissertation, the Einstein 
convention 

N 

a jbj = (2*4) 

j=i 

will be used to simplify the form of the equations. Here N is the number of physical 
dimensions in the problem. 
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Rewriting the nonlinear term in equation 2.3 using the identity 

, (UfrUfc)^ 

UkUj'k = CjklUlUlk + 




yields 


u j,j — 0 
,p , u k u k , 
P 


Ujj + tj k iuiuj k + (“ + 2 — vu i' kk 


(2.5) 

(2.6) 

(2.7) 

( 2 . 8 ) 


where lJj is the vorticity component in the Xj direction, and t 3 ki is the alternating 
unit tensor 

(1, if (j,k,l) = (1,2,3), (3,1,2), or (2,3,1) 

-1, if (j,M) = (3,2,1), (1,3,2), or (2,1,3) (2.9) 

0, if j = k, j = /, or k = /. 


tjkl 




Taking the curl of equation 2.8 and using the identity 


ejki<t>,ik = 0 


(2.10) 


where <f> is any scalar, yields the vorticity form of the Navier-Stokes equations. 


T ejfc(Cimn( lt n tL, m) 1 *: — vl -^j,kk 


( 2 . 11 ) 


U)j = tjklUl' k . 


( 2 . 12 ) 


These equations will be solved along with the set of passive scalar equations 


Cj, t + UkCj'k — K G Jy 


kk 


(2.13) 


where C 3 is a scalar concentration and n is the scalar diffusion coefficient. Since 
equation 2.13 is linear and homogeneous in Cj the magnitude of Cj is arbitrary, 
hence normalization is irrelevant. In addition, because of the form of equation 2.13 
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Cj may be offset by any arbitrary constant. This affords wide latitude in the choice 
of initial conditions for the scalar field. 


2.3.2 Stretching Grid 

In order to maximize the available resolution for a given number of grid points while 
at the same time keeping the vorticity at the edges of the resolved box small to satisfy 
the asymptotic matching condition, a growing uniform grid is used in the non-periodic 
x 2 direction. To implement this, the coordinates in each direction are rescaled to an 
interval of length 27r, the “natural” interval for Fast Fourier Transforms. 


6 = 


27 ra-j 

~17 


6 (<) = 


2irx2 

m 



(2.14) 


where L\ and L 3 are the fixed box lengths in the periodic Xi and X 3 directions re- 
spectively, and L 2 (t) = 12(0 “ Vi(0 is the time varying box size in the aperiodic x 2 
direction. For convenience it will be assumed that Y 2 (t) = —Y^(t). 

Applying these coordinate transforms to equations 2.11, 2.12, and 2.13 yields 


Wi. 


C,t- 


l 2 


d~ €jkl {€lnq'UqUJ n ) >m ^mk — A/ m A m fc 

(2.15) 

— ^jkl^l,m^mk 

(2.16) 

Cj, 2 + ‘U'kC'j,m ^mk = j t kl^lm^mk 

(2.17) 

_ f 2 t v/Lj, if j=k; 

(2.18) 

( 0 , otherwise. 


The second term on the left hand side of equation 2.15 and the similar term on 
the left hand side of equation 2.17 are due to grid stretching. They may be absorbed 
into the nonlinear terms by defining a grid-relative velocity u* = uj — L 2t t(t)^ 2 6 2 j/27r. 
In terms of this modified velocity the equations of motion become 


&j,t + A m * + ')u}jL 2 ,t! L 2 — vu h kiAi m A m k 


(2.19) 
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Wj — — ^jkl u l y m^mk ( 2 . 20 ) 

Cjj “1“ 'U'i c Cj t rn^km = &Cj % kl Aim A m fe (2.21) 

A . k = \ 2n/Lj ' l/ J=k; (2.22) 

( 0, otherwise. 

7 = ( 1, j =1,3; (2.23) 

( 0, otherwise. 


The remaining grid stretching terms in equation 2.19 are in a form that may be readily 
Fourier transformed. 

In domains I and III, the terms in equation 2.19 are identically zero, and the 
left hand side of equation 2.20 goes to zero. Equation 2.20 will be solved to ob- 
tain asymptotic solutions for the velocity in domains I and III and total velocity in 
domain II. 


2.4 Numerical Method 

2.4.1 Transformed Equations 

The governing equations 2.19 through 2.21 are solved in domain II using a standard 
Fourier pseudo-spectral technique treating the vorticity as periodic in the aperiodic 
direction £ 2 - Note that all of the terms in equation 2.19 can be treated as periodic 
since by construction — » 0 at the top and bottom of the computational box. The 

same holds true for the scalar convection equation 2.21. The governing equations are 
advanced in wave space with nonlinear terms which are calculated in physical space 
at each timestep. 

Equations 2.19, 2.20, and 2.21 Fourier transformed in all three directions become 
Uj, t + (^2i + isu J k k kiAi m \ mk ) + ik m ^{e jkl (ti ng u*u n )\ km ) = 0 (2.24) 

l 2 


LOj — ICjjfcjU, k m A mk — ItjklUlkmAmk 

Cj lt T K.Cjk k kiAi m A mk d" F(Cj' k u m A mk ) 0 


(2.25) 

(2.26) 
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f 2x/ Lj , if j=k; 

[ 0, otherwise. 

_ 1 1, if j=l,3; 

( 0, otherwise. 


(2.27) 

(2.28) 


where kj is the wavenumber in the direction, hatted quantities are Fourier trans- 
forms of the corresponding physical vector, and T{) is the Fourier transform operator. 

The diffusion term in equation 2.26 as well as both the diffusion term and the grid 
stretching term which appear in equation 2.24 are absorbed into the time derivative 
terms by use of integrating factors, yielding the basic set of equations in wavespace. 


1 

W) 


(G(t)uJj') j 4“ ikmIE(£jkl(€l r iq‘UqUJ n ')A rn k') 0 


(2.29) 


LmAmk — ^^■jkl’^l^mAmk 


1 


H(t) 


(H(t)C j ) >t + f(ulA mk C hm ) = 0 


G(t) = Ll(t)exp(vkik k f A /m A mfc dr) 

Jt o 


H(t) = exp (Kk t k k f Ai m A m kdT) 

Jt 0 


hjk = 


f 27 r/Lj, if j=k; 

( 0, otherwise. 




1, if j=l,3; 
0, otherwise. 


The solution procedure for this set of equations is as follows. 


(2.30) 

(2.31) 

(2.32) 

(2.33) 

(2.34) 

(2.35) 


• The periodic part of the velocity field is calculated from the periodic vorticity 
in wave space for domain II. 


• The velocity and vorticity fields are transformed into physical space in the 
aperiodic direction. 
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• An aperiodic field is added to the periodic velocity field to form the total, 
aperiodic, velocity. 

• The velocity and vorticity fields are transformed into physical space in the 
remaining two directions and the nonlinear terms are formed. 

• The nonlinear terms, which are periodic by virtue of the vorticity going to zero 
at the edges of the box, are then transformed in all directions back into wave 
space. 

• The nonlinear terms are dealiased. 

• The nonlinear terms are used to advance the vorticity field in time. 

2.4.2 Asymptotic Matching of Velocities 

The crux of the procedure outlined above is the calculation of the aperiodic component 
of the velocity field. This is calculated from the curl of 2.20 which is 

u j,kl^lm^mk = ^jkl^l.m^mk- (2.36) 

Only one component of Uj need be solved for, since the remaining two components 
can be constructed from continuity and the definition of vorticity. It is convenient to 
solve for the cross-stream velocity u 2 . The equation for u 2 is transformed in the two 
periodic and £3) directions, but not in the non-periodic (£ 2 ) direction giving the 
second order ordinary differential equation 

-^2^2,22 + (“jF| + )^2 = -J^h&klt2lm&m- (2.37) 

For K 2 = Ll(kl/Ll + k\j L\) > 0, the general solution to equation 2.37 is 

u 2 = + Aexp (-M^- + |f )*&) + Bexp(-L 2 (-^ + ^|)^ 2 ) 


(2.38) 
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where u % is the particular solution of equation 2.37 (zero in domains I and III where 
the vorticity is zero, calculated numerically from the vorticity field in domain II), and 
A and B are functions of K which are determined by asymptotic matching of the 
velocity and velocity gradient. 

The solution for the velocity must be bounded at £2 — ± 00 , hence in the three 
domains 

Aje h domain I 
62 = < uf + Ajje h + Bue~ K( ‘ '- 2 domain II (2.39) 

, Biiie~ K( - 2 domain III. 

Matching the velocities and velocity gradients at the domain boundaries £2 = ±7r 
yields four algebraic equations in four unknowns 

Aje-* K = fif(-w) + Aut~ vK + B u ^ k (2.40) 

KAit~* K = uj 2 (- tt) + KAne~* K - KB u e* K (2.41) 

Bjue~ wK — uf(7r) + Aije* K + Z?//e _,rK (2-42) 

- KB lu e ~ tK = fij 2 ( x) + KA n e* K - KB u t~ vK . (2.43) 

Solving for Ajj and B// gives 


Ma ']^ + «r( 7r ))e^P(^6)) 

(2.44) 

(2.45) 

Solutions for u 1 and U 3 are constructed from the incompressibility condition and the 
definition of vorticity. 

If K = 0, equation 2.37 is replaced by the two-direction transform of 2.20 with 
ki = k 3 = 0 



&j = tj 2 lUl,2^22- 


(2.46) 
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Integrating in the £ 2 direction yields 

(2.47) 

where D 3 are constants which are determined by the mean velocity in each direction. 

Solutions for u \ and U 3 are constructed from the resulting solution for u 2 > n do- 
main II for each (ki, k 3 ) pair and the total solution for the velocity is transformed to 
physical space in the periodic directions. The cross product of the velocity and the 
vorticity is taken to form the periodic nonlinear term. The nonlinear term is then 
transformed back into wave space to advance the vorticity and scalar fields in time. 



2.4.3 Time Advance 

The time advancement method used for all dependent variables is the second order 
Runge-Kutta scheme 


& 

II 

e- 

(2.48) 

$" +1 = $ n + A</($ n ;t n ) 

(2.49) 

$" +1 = <J>" + ^(/($";<») + /($” +1 ;< n+1 )) 

(2.50) 


where n is the time index. 

This particular scheme was chosen because it minimizes memory requirements for 
a given simulation size. It requires only two copies of the three vorticities and each 
scalar (one field at t = t n and one at t = t" +l ) to be stored at any given time. It also 
allows for the use of the random-shift dealiasing scheme described in section 2.4.5 and 
appendix C, which helps make maximum use of the available computational power. 

Applying this scheme to the governing equations 2.29 and 2.31 yields 
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C " +1 


j£rCj + f (jg^KA mk C hm ) 


(2.52) 


The precise forms of G n and H n depend on the way in which the grid is stretched. 
In the present code, grid stretching was taken to be piecew'ise linear, hence 


(2.53) 



and 


H n 

tfn+l 


/ , k\ H i lx H \ j\\ 

= exp(-47T k{ At + — -(— - + j2 Ai ))• 


This is the form of the equations implemented. 


(2.54) 


2.4.4 Accuracy and Stability 

As implemented, the code has standard Fourier spectral accuracy in all three spa- 
tial directions and is second order accurate in time. The second order Runge-Kutta 
timestepping scheme used is “weakly” unstable (it is unstable in the absence of vis- 
cosity, though only mildly so for small At). However in the presence of even a small 
viscosity the method is stable for a range of At. 

This can be quantified by applying the method to the one dimensional linear 
convection-diffusion equation 


U,t -I- CU' X = WU 'lx 


for a single Fourier component 


u 


= C(t)exp( 


2 nikx^ 

TV Ax ' 


(2.55) 


(2.56) 
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Figure 2.3: Neutral stability curve for second order Runge-Kutta timestepping algo- 
rithm applied to linear con-vec-tion-diffusion equation u^ t Acu iX = i su iXX . Coordinates 


are B = 


27 ri/k 
cN Ax 


and r = 


2rkcAt 
N Ax 


Substituting equation 2.56 into equation 2.55 yields 


/ 27 rikc 47 v 2 k 2 i>\ 

~ \iVAx + iV 2 A x 2 J 

= Lu. 


(2.57) 


Applying the timestepping algorithm from equations 2.48 through 2.50 to equation 
2.57 gives 

u n+1 = (l + LAt + ^-\ u n . (2.58) 


The timestepping method is considered stable if u n remains finite as n — > oo. This 
will be true if 


u" +1 

2 

u n 



= 1 + LAt + 


I 2 At 2 


(2.59) 


It is convenient to define the dimensionless variable 


_ 2xkcAt 2irk 
r = ~NAx~ = ~N~ 


(2.60) 
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and the dimensionless parameter 

D _ 2irvk 

= cTVAx 

where CFL is the Courant-Friedrichs-Lewy number. In terms of these dimensionless 
quantities 

LAt = -{B + i)r (2.62) 

and the stability condition from equation 2.59 expands to 

( g2 + Hl r 3 - B{B 2 + \ )r 2 + 2 B 2 r - 2B < 0. (2.63) 

The neutral curve for equation 2.63 is shown in figure 2.3. The equation is third 
order in r (hence third order in At) and no simple analytic stability criterion exists. 
In order to choose a A t, therefore, a simple Newton method solver was used to solve 
equation 2.63 numerically for r(B ) at each timestep. Note that for a given timestep 
the parameter B is a quantity with known limits (at k = 1 and k = N/2). At each 
timestep r(B) was evaluated twice, once for k = 1 and once for k = N/2. The 
minimum of the two resulting values for A f, multiplied by a factor of C = 0.8, was 
used for the next timestep size. 

The factor C < 1.0 was necessary to ensure time accuracy in addition to stability 
of the method. Several test runs with identical initial conditions (the same as the test 
case discussed in section 2.8.3) but with C varying from 0.4 to 1.0 were performed. 
The value C = 0.8 was found to be sufficiently small for the solution to converge to 
a solution independent of C . 


(2.61) 


2.4.5 Alias Control 


Following the approach of Rogallo [34], aliasing was controlled with a combination 
of high wavenumber masking and phase shifting. All modes with wavenumbers such 


that 


kl k\ k\ 2 

— — H — H — > - 

N 2 ^ N 2 T N 2 g 


(2.64) 
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are set to zero. This eliminates all two- and three-dimensional aliasing, leaving only 
the one-dimensional aliasing term in each direction. The one-dimensional alias errors 
are dealt with by phase shifting the data a random fraction of a grid cell width in 
each direction at every other time substep followed by a further shift of exactly one 
half a grid cell width in each direction at the subsequent time substep. This random 
shifting cancels the aliasing error to second order in time, the same order as the 
time advancement algorithm. Using this dealiasing technique instead of a perfect 
(2/3 rule or multiple phase shift) technique reduces both memory requirements and 
operations per timestep. See Appendix C for a more complete discussion of aliasing 
and dealiasing for Discrete Fourier Transforms. 


2.5 Code Implementation and Data Management 

The code described above has been implemented on an Intel iPS< 7860 supercomputer. 
Also known as the “Hypercube”, the iPSC/860 used is a massively parallel computer 
which has 128 computational nodes each consisting of a 40 megahertz Intel i860 
processor with 8 megabytes of random access memory (RAM). The code is written 
in VECTORAL, a high level programming language developed by Wray [57] w r hich 
facilitates the handling of large data sets. 

The architecture of the machine dictated the implementation of the code. The 
whole simulation, executable and data, had to fit in the distributed RAM on the 
machine. This requires careful division of the data into useful subsets which reside 
separately on each processing node. The method described is natural to implement 
on a machine with such an architecture since nearly all of the operations at any given 
point in the algorithm are performed independently on separate parts of the data set. 
This leads to a simple division of tasks among the multiple processors, making for a 
fast and efficient code. 

In order to perform a transform on a line of data, a given processor must have 
all the data in the transform direction for that line. To achieve this, the data was 
structured for a two pass method, with each of the N computational nodes storing 
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Table 2.1: Code Structure 

Pass 1 

1 . Transform nonlinear terms into wave £2 space, shifted mesh. 

2 . Unshift nonlinear terms. 

3 . Advance governing equations for substep. 

4 . Save data and/or stop if necessary. 

5 . If first R-K substep, calculate At, dealiasing shifts, integrating factors, etc. for 
timestep. 

6 . Zero “oddball” ( kj = Nj/2) wavenumber contributions. 

7 . Zero high wavenumber modes to eliminate multi-dimensional aliasing. 

8. Shift data for one-dimensional alias control. 

9 . Form needed £2 derivatives of data. 

10 . Transform data into physical £2 space, shifted mesh. 

11. Execute Pass 2 . 

Pass 2 

1. Transpose data to (£1 , £3) planes. 

2. Calculate additional velocity and form corrected velocity terms. 

3 . Form needed £i and £3 derivatives of data. 

4 . Transform data into physical £1 and £3 space, shifted mesh. 

5 . If second R-K substep, calculate max (^) for CFL stability requirements. 

6. Form nonlinear terms in physical space, shifted mesh. 

7 . Transform data into wave £1 and £3 space, shifted mesh. 

8. Transpose data to (£1 , £2) planes. 


9 . Execute Pass 1. 
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and manipulating l/N of the total data set during each pass, and swapping data 
between passes. 

For the first pass, each computational node holds all the 6 and 6 data for 1/A r 
of the £3 planes. All £ 2 derivatives and transforms are evaluated during this pass. For 
the second pass, the data is swapped between computational nodes so that each node 
now holds all of the £1 and £3 data for l/N of the £2 planes. All £1 and £3 derivatives 
and transforms are evaluated during this second pass. The data is then swapped back 
to its original configuration and the governing equations are advanced in time. Table 
2.1 lists the operations executed during each time substep. 


2.6 Boundary Conditions 

Periodic boundary conditions are imposed in the streamwise £1 and spanwise £3 di- 
rections 

$(£1 + 2 ?r, £2, £3; ^ ) = $(£1,6,6; t) ( 2 - 65 ) 

$(6 , 6, 6 + 2tt; t ) = $(6 ,6,6 ; i) (2.66) 

where $ is any dependent variable. 

The boundary conditions imposed on the vorticity field in the cross-stream direc- 
tion £2 are that the vorticity is zero outside the resolved box 

u>j = 0 for 6 < 0 or 6 > 27 t ( 2 . 67 ) 

and that the velocity perturbations go to zero at infinity 

1 - ( 2 . 68 ) 

The boundary conditions imposed on the scalar field in the cross-stream direction 
are that each scalar concentration be the same at the top and bottom of the resolved 


box 


Q(£i,0,6) = Q(6,2t,£ 3 ). 


( 2 . 69 ) 
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This is most easily satisfied by picking initial scalar concentration distributions that 
go to zero in the freestream. 


2.7 Initial Conditions 

2.7.1 Vorticity Field 


The initial conditions for the vorticity field of the time developing wake consist of a 
Gaussian mean streamwise velocity profile 

Ui = 1 — Au(£ 0 e~ ClX 2 

(2.70) 

u 2 = u 3 = 0 

(2.71) 

which gives the mean vorticity profile 


ijJl — Us'2 ~ 0. 

(2.72) 


= -2 Cl Au^ 2 e~ c ^ (2.73) 

The centerline velocity defect, Au| 0 , was chosen to be 0.692 and Cj was chosen to 
be 0.69315. This particular profile was used in the experiments of Sato & Kuriki [39] 
and Corke, Krull, &: Ghassemi [14], and in the computations of Chen, Cantwell, & 
Mansour [10]. It gives an initial wake halfwidth of 6o = 1.0 and an initial Reynolds 
number based on halfwidth of Re), = 0.692/i/. 

For this mean velocity profile, the relationship between the various possible wake 
width length scales is as follows: 

b 0 = 1 (2.74) 

f>o — J (1 - u\)dx 2 = 2 Au<£ 0 6o = 1.473&0 (2.75) 

6eo = J°° «i(l - ui)dx 2 = ^ f 1 - Au^obo = 0.75246 o . (2.76) 
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Though the wake halfwidth, bo , is used here as the reference scale (see figure 2.2 on 
page 22), the other scales given above will be used later when comparing results to 
da: • from other studies. 

Small disturbances which were periodic in the streamwise and spanwise directions 
were added to the mean flow. 


ft(6, 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


6 , 6;0 ) = ft 

Real[e 100^100 exp(mfi ) 

A _/ Q fi + <^010 ^ 

£oio‘*oioexp(t — ) 

A -i ■ a 6 + \ 

£ooi“ooi exp(z - ) 

. 100 A 100 / 


£ 100 fl 100 (exp(i(a^i - fiiz + <£ 100 )) - exp(t'(a£i + /?6 + ^ 100 ))) 

,.«6 +^3 + r u 

) - exp(t )) 


£ 010 fl 010 (exp(r^— + + ^ 01 ° 

e M 1 n 001 ( exp (« a ^~'^ 3 + ^° 01 - 


,.a6 + 06 + *°°\* 

) ~ exp(z )) 


...] 


(2.77) 


where ft is any of the three vorticities, ft = ft(6) is the mean flow, and ft = ft(^) is a 
disturbance eigenfunction determined from linear theory as described in appendix B. 
The quantities o and 0 are the streamwise and spanwise wavenumbers of the funda- 
mental mode respectively, and <j> is a disturbance phase angle. Subscripts indicate a 
two-dimensional disturbance, superscripts indicate a three-dimensional disturbance, 
and the positions of the Ts indicate the streamwise wavelength of the mode. As an 
example, £ 001 is the disturbance amplitude of the three-dimensional disturbance at 
the sub-subharmonic streamwise wavelength. 

The disturbance phases, <f > , are all defined with respect to the two-dimensional 
fundamental disturbance, with <j> = 2tt corresponding to a physical shift of one fun- 
damental wavelength. For the two-dimensional disturbances, the phase is measured 
from the first zero (£i = 0) of the fundamental velocity disturbance to the first zero 
of the longer wavelength velocity disturbance as shown in figure 2.4. Figures 2.5 
and 2.6 show the phase for the three-dimensional fundamental and three-dimensional 
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subharmonic disturbances respectively. In each case, the phase is measured from the 
first maximum of the two-dimensional vorticity fundamental disturbance to the first 
crossing of the maxima of the pair of oblique vorticity disturbances. 

Hence <^ 010 = 0 represents the phase of a three-dimensional subharmonic of the 
vorticity with its first crossing points aligned with the maximum of the two-dimension- 
al fundamental. Similarly, d> 010 = 7r/4 represents a three-dimensional subharmonic 
with its crossing points aligned with zeros of the two-dimensional fundamental. 

All of the flows simulated had two-dimensional disturbances at a wavelength equal 
to the fundamental (Karman) mode superimposed on the mean (eioo > 0). The three 
dimensional disturbances were pairs of oblique waves with equal and opposite span- 
wise wavenumbers oriented 60° to the spanwise direction (f 3ja — tan(60°)). This 
choice of three-dimensional modes was motivated by the stability analysis of Flem- 
ming [16] which suggested that the most unstable three dimensional modes should 
be wave pairs at angles near 60° at the subharmonic wavelength. The experimental 
observations of Corke, Krull, & Ghassemi [14] supported this analysis. 

The amplitude for each disturbance eigenfunction was chosen such that the inte- 
gral over x 2 of the magnitude of the disturbance velocity eigenfunction for the given 
mode was equal to 0.02£/o&o- This initial magnitude was found to be small enough 
for the initial wake growth to be within the linear regime. At the same time it was 
large enough to allow the wake to enter the non-linear growth regime without undo 
expenditure of of computational time. 
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Figure 2.5: Disturbance phase for three-dimensional fundamental. 

2.7.2 Passive Scalar Field 

Though the method allows for carrying an arbitrary number of passive scalars, in 
practice only one scalar was carried in the simulations. The initial passive scalar 
concentration at each point was taken to be the magnitude of the vorticity (the 
enstrophy density) 

Cl = (u> k 0Jk)> ■ (2.78) 

Since the vorticity perturbations were small, the initial scalar field is very nearly the 
square of the mean vorticity field 

Cl « 2c 1 Au< E oMe- Cl < (2.79) 

After the simulation was initiated, the vorticity and scalar fields were allowed to 
evolve independently. 
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Figure 2.6: Disturbance phase for three-dimensional subharmonic. 


2.7.3 Grid Stretching Rate 

The grid stretching rate, L 2 ,t, was initially set to follow a L 2 ~ {t — to ) 1 ^ 2 growth 
curve to match the expected long-time growth of the wake. This gives 


I 2 ,t(0 


L\{t = to) 

2toL2{t) 


(2.80) 


where t 0 is a virtual time origin which depends on the initial Reynolds number. 

The growth rate was periodically adjusted manually as the simulations progressed 
to keep the magnitude of the mean vorticity at the boundaries of the resolved box 
below a fixed percentage of the maximum mean vorticity magnitude. This manual 
adjustment was primarily necessary during the initial stages in the development of 
the wakes when the growth was significantly different than the asymptotic t 1 ? 2 curve. 
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2.7.4 Flow Reference Convention 

In order to facilitate references to the various wake flows discussed throughout the 
remainder of this document, it is useful to introduce a standard naming convention 
which will present all the relevant information about a given flow in a compact form. 
From this point forward, the simulations presented here will be referenced with a tag 
of the form 

( 2 . 81 ) 

where ( R ) is the initial wake halfwidth Reynolds number, Re b , (A) is the angle of 
the three-dimensional disturbances with respect to the spanwise direction (zero if 
the flow is two-dimensional), and (2D) and (3D) are three-element strings which 
indicate which disturbance modes are present and gives their phases with respect to 
the two-dimensional fundamental disturbance. 

The possible elements in the strings (2D) and (3D) are x, which indicates the 
given mode is not present, and 0, tt/ 4, or 7t/2 which indicate that the given mode is 
present and is at the corresponding phase with respect to the fundamental. The loca- 
tion of the element in the string gives the wavelength of the disturbance as outlined 
above for e and 4>. The first slot indicates a disturbance at the fundamental stream- 
wise wavelength, the second indicates a disturbance at the subharmonic streamwise 
wavelength, and the third indicates a disturbance at the sub-subharmonic streamwise 
wavelength. 

As an example, the tag 346(60)o|* refers to a wake with a Reynolds number of 
Re b = 346, and disturbances corresponding to a two-dimensional fundamental (which 
must be at a phase of 0 to itself), a two-dimensional subharmonic with <^ 0 io = * /2, 
and a pair of 60° oblique disturbances with <t> 010 = tt/4. 

In addition, the wildcard character ‘(?)’ will be used to indicate that the given 
parameter can take on any appropriate value. For example, the tag 346(60) o i]. 
refers to all Re b = 346 wakes with a two-dimensional fundamental, a two-dimension- 
al subharmonic with <£ 0 io = w/4, a P air of 60 ° oblige disturbances at any of 
the calculated phases. The tag (?)(60)££ refers to all wakes with a two-dimensional 
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fundamental, a pair of 60° oblique disturbances with <f> 010 = 0, and any of the initial 
Reynolds numbers simulated. 


2.8 Code Validation 

The code has been validated with three sets of test cases. The first two tests, pure 
diffusion and linear disturbance cases, compare results from the code to linear flow 
solutions. The third is a comparison of results for identical initial conditions between 
the present code and the well tested and generally accepted code of Spalart, Moser, 
& Rogers [44] The results of these test cases are outlined below. 

2.8.1 One Dimensional Diffusion 

For two-dimensional parallel flow in the aq direction, the Navier-Stokes equations 
reduce to the one-dimensional diffusion equation 


« 1 ,« = ^ 1 , 22 - 


(2.82) 


This equation admits analytic solutions for the diffusing parallel wake 

“■ “ 1 - ( 2 - 83 ) 

and for the diffusing parallel mixing layer 


uj — 2erf( 


Xi 


(4 u(t - f 0 )) J 


r) - L 


(2.84) 


To check agreement with the analytic solutions, four two-dimensional test cases 
were run: two parallel wakes started from a Gaussian velocity profile, one with grid 
stretching turned off and one with grid stretching turned on; and two parallel mixing 
layers started from an error function profile, again with grid stretching off in the first 
and on in the second. Each case was run on a 4 x 128 x 4 grid on two computational 
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Figure 2.7: Halfwidth vs. time for Re b = 35 diffusing parallel wake. L 2 ( 0) = 10.0. 
□ : Computed wake on fixed grid, l 2 ,t(0 = 0.0. O: Computed wake on growing grid, 
L 2 ,t{t) = 0.125. : Analytic solution. 



Figure 2.8: Vorticity thickness vs. time for Re e = 50 diffusing parallel mixing layer. 
I 2 (0) = 10.0. □: Computed mixing layer on fixed grid, L 2yt (t) = 0.0. O: Computed 
mixing layer on growing grid, L 2 j(t) = 0.18. : Analytic solution. 
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Figure 2.9: Centerline velocity defect and total x 2 momentum vs. time for Ret, = 35 
diffusing parallel wake. L 2 (0) = 10.0. □: Computed wake on fixed grid, L 2it (t) = 0.0. 

O: Computed wake on growing grid, L 2 i t (f) = 0.125. : Analytic solution for 

centerline velocity defect. : Analytic solution for total x 2 momentum. 

nodes. Two computational nodes were not required by the problem size, but were 
used to verify internode message passing. 

The numerical results and the analytic solutions are plotted in figures 2.7, 2.8, 
and 2.9. The wake flow halfwidth, 6, is defined as half the width at half the centerline 
defect velocity. The mixing layer vorticity thickness is defined as the freestream 
velocity difference divided by the maximum mean vorticity (in this case the vorticity 
at the centerline). 

The agreement between the numerical results and the analytic predictions for 
t < 50.0 is excellent for both wake cases and for t < 225.0 for both mixing layer 
cases. Past a time of t = 50.0 for the wake and t = 225.0 for the mixing layer the 
simulations with fixed grids begin to deviate significantly from the analytic solutions 
as the vorticity field outgrows the resolved domain. This violates, benignly, the 
requirements of the numerical scheme. 

The excellent agreement for even very late times for the cases with the grid growing 
to keep pace with the spreading vorticity indicates that the diffusion portion of the 
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code, including the implementation of the growing grid, works properly. Over the 
course of the runs the wake simulation grid grew by a factor of nearly five and the 
mixing layer grid grew by a factor of over nineteen. These are rather stringent tests 
of the diffusion part of the implementation. 

2.8.2 Linear Growth Rate 

The second verification test for the code was a comparison between eigenmode growth 
rates for a computed solution of the Navier-Stokes equations using the presem code 
and the eigenmode growth rates predicted by linear stability theory. The numerical 
simulation started from a Gaussian mean profile perturbed by a very low amplitude 
most-unstable eigenfunction (as predicted from linear theory). The eigenmode growth 
rate was calculated by assuming that each eigenmode can be represented by a wave 
having complex growth rate ac and complex modeshape f/j. In particular, for the u\ 
velocity component 

u 1 (a,x 2 -,t) = U 1 (a,x 2 )e-' act . (2.85) 

Taking the logarithm of equation 2.85 and solving for the real part of the result yields 
an equation for the growth of the disturbance 

acit = Real[ In tti(a, x-i ; t) — In U\(a, X 2 )]. (2.86) 

Hence, the growth rate, ac,, is the rate of change of the real part of of the logarithm 
of the given dependent variable. 

Figure 2.10 shows the calculated linear growth rate and the predicted linear growth 
rate for the test computation. The numerical simulation was performed on a 128 x 
128 x 4 grid on 4 processor nodes. The flow Reynolds number was set at Re t = 1384 
to reduce the effect of viscosity on the results as much as possible without requiring 
excessive resolution. Agreement between prediction and the computation is excellent. 
Deviations in growth rate near t = 0 are due to small errors in the disturbance 
eigenfunction incurred while mapping the disturbance from the nonuniform numerical 
mesh used by the linear stability solver code to the uniform numerical mesh used in 
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Figure 2.10: Magnitude of two-dimensional fundamental disturbance vs. time for 
parallel wake. : Linear stability prediction. O: Computed flow. 

the simulation code. The slight dropoff in the growth rate at late times for the 
computed flow is due to viscous spreading of the mean flow, which is not accounted 
for in the linear theory. 

2.8.3 Comparison with Accepted Jacobi Polynomial Code 

The final test consisted of a comparison between a numerical simulation using the 
present code and a numerical simulation started from identical initial conditions using 
the accepted code developed by Spalart, Moser, &: Rogers [44] (hereafter referred to 
as the SMR code). 

Figure 2.11a shows the three-dimensional vorticity magnitude for a Reb = 346 
wake started from a Gaussian mean velocity profile with a two-dimensional distur- 
bance at the most unstable (Karman) wavelength and a three-dimensional disturbance 
at the subharmonic wavelength and sixty degrees from the spanwise direction. For 
this case grid grow r th was active with Z/ 2 ,t = 0.1, which translates to a 50% increase 
in grid size from the initial grid to the grid at the time shown. This condition should 
provide the most stringent test of code performance. 
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(a) 



(b) 


Figure 2.11: Enstrophy Density isosurface from direct numerical simulations of 

346(60)g°* wake. |w| = 0.7. (a) Present code, t = 103.3. (b) SMR code, t = 104.4. 
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Figure 2.12: Wake halfwidth vs. time for direct numerical simulation of 346(60)£°£ 
wakes. : Present code. □ SMR code. 

Figure 2.11b shows the identical flow calculated using the Rogers and Moser code. 
Some differences between the two calculations are expected due to differences in the 
computational grid and time advance algorithms. Despite this, the magnitude and 
shape of the structures in both codes are very similar even at the advanced time shown. 
Such excellent agreement is a strong indicator that the code functions properly. 

As a check of the resolution of the test simulation, two variants of the Taylor mi- 
croscale, A t (calculated from all three velocity components) and A71 (the traditional 
microscale calculated from only the streamwise velocity) were computed as func- 
tions of £ 2 (see appendix A for definitions of the microscales used). Both measures 
gave a minimum nondimensional microscale length of approximately 0.9 (Xj = 0.83, 
Aji = 0.87) at the time shown in figure 2.11a. That length corresponds to a minimum 
of approximately 5 grid points in any direction, thus it is reasonable to conclude that 
the flow is sufficiently resolved. 

Figure 2.12 is a plot of the square of the flow halfwidth versus time for the same 
two simulations. This represents a good measure of how well the codes match mean 
flow characteristics. There is no significant difference between the results obtained 
from the two codes. The small deviations at late time are due to small difference in 
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flow time and the fact that the vorticity field is beginning to outgrow the well resolved 

portion of the fixed mesh in the SMR code. 

This test also gives a rough comparison of performance between the two codes. 
The iPSC/860 run required 5 hours on 16 of the available 128 nodes to run to the 
time shown. The Cray simulation took approximately 15 cpu hours spread out over 
2 weeks of run time on a single processor of an 8 processor Cray-YMP. 

2.8.4 Behavior Near the Matching Boundaries 

The fundamental assumption of this method is that the vorticity remains confined 
to the resolved domain (it is compact) and therefore the vorticity magnitude at the 
matching boundaries remains negligible. It is important to know if this assumption 
is valid in cases of interest, and what consequences arise when these assumptions are 
violated. 

Figures 2.13 and 2.14 show the behavior of the vorticity magnitude, normalized 
by the maximum vorticity magnitude in the flow, near the matching boundaries for 
two wake simulations. Figure 2.13 is data from the low Reynolds number two-dimen- 
sional diffusing w r ake discussed in section 2.8.1. Figure 2.14 is data from the higher 
Reynolds number strongly three-dimensional wake discussed in section 2.8.3. The 
sample line for the data in this figure was chosen to be the one on which the high 
enstrophy density regions of the wake pass closest to the edge of the box at late times. 
This line of data represents the worst case condition. 

At early times in each flow, the behavior of the vorticity magnitude near the 
matching boundaries fit very well with the assumptions of the numerical method 
(max \i0,\ edge < max \u)j\jiow)- Even at late times the flow near the boundary remains 
well behaved. The small amplitude ripples near the matching boundary for the late 
time of the Re b = 346 wake are due to the slight mismatch in the vorticity magnitude 
between the top and bottom edges of the resolved domain. Attempting to resolve this 
small jump using a finite set of basis functions results in Gibbs phenomenon waves 
near the jump point. 



CHAPTER 2. NUMERICAL METHODOLOGY 


50 



Cross-stream grid point 


Figure 2.13: Behavior of relative vorticity magnitude near the matching bound- 

aries for parallel diffusing wake, Rei, = 35. □: t = 1.9. O: t = 168.9. Matching 
boundary occurs between grid point 0 and 1 (left edge of plot). 
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Figure 2.14: Behavior of relative vorticity magnitude near the matching bound- 

aries for three-dimensional wake, Re = 346. □: t = 2.8. O: < = 103.3. Matching 
boundary occurs between grid point 0 and 1. 
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These small ripples are not a significant source of error in the simulations for a 
number of reasons. First, the magnitude of the Gibbs phenomena ripples are pro- 
portional to the magnitude of the vorticity mismatch across the edges of the resolved 
box. Thus so long as the vorticity mismatch is kept small, the magnitude of the 
ripples will be small. Second, the Gibbs phenomena ripples represent both positive 
and negative deviations in the vorticity field (note that figure 2.14 shows magnitude). 
Since the ripples are confined primarily to the edges of the resolved box, far from 
the vorticity containing region of the wake, the integrated (Biot-Savart) effect of the 
ripples is much smaller than the magnitude of the ripples would suggest. Finally, 
the Gibbs phenomena ripples occur at the highest wavenumbers, which are strongly 
damped in a well resolved viscous flow. Here the magnitude of the Gibbs waves are 
less than 0.2% of the maximum mean vorticity at their worst (cut through widest part 
of the wake). This magnitude is easily small enough to keep from having a significant 
influence on the flow. 



Chapter 3 

Two-dimensional Simulations 


3.1 Motivation 

In order to have a reference with which to compare the three-dimensional wake calcu- 
lations presented in the next chapter (starting on page 93), a set of two-dimensional 
plane wakes was simulated with a variety of initial conditions and Reynolds numbers. 
The effect of these parameters on the evolution of two-dimensional structures in a 
wake and on the mean velocity profile is examined in this section. 

3.2 Simulation Parameters 

Two-dimensional simulations require a relatively small investment of computational 
resources, and are therefore ideal for parametric investigations. Simulations with a 
variety of combinations of disturbance wavelength and phases, and flow Reynolds 
numbers were run. All of the simulations were initiated from a Gaussian mean wake 
profile upon which a most unstable anti-symmetric (Karman) disturbance is superim- 
posed. Combinations of subharmonic and sub-subharmonic disturbances were added 
to the initial fields in some of the simulations. Wake Reynolds numbers based on 
initial wake halfwidth and centerline defect velocity were varied from Re*, = 69 to 
Reb = 2768. The initial passive scalar field in all the simulations was set equal to the 
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Figure 3.1: Two-dimensional disturbance phasing 
initial enstrophy density at each point 

Cj = (u>k^k) 2 ■ ( 3 - 1 ) 

and the Prandtl number 

Pr = - (3-2) 

K 

was set to a value of 0.7. Table H.l in appendix H lists all the two-dimensional wakes 
simulated and their parameters. 

Figure 3.1 gives a graphic summary of the disturbance wavelengths and the mean- 
ing of the disturbance phases used in the two-dimensional simulations. Initial distur- 
bances were combinations of most-unstable eigenmodes at the fundamental, subhar- 
monic, and sub-subharmonic wavelengths. Note that the phase is always in reference 
to the wavelength of the fundamental, so (for example) a phase of 7t/2 results in a 
shift of one quarter of a fundamental wavelength. See section 2.7 for a more general 
description of the disturbance functions. 

3.3 Evolution of the Two Dimensional Wake 

3.3.1 Initial Development 

For all of the two-dimensional wakes studied, the early evolution followed the same 
general pattern. This is typified by figure 3.2, which shows contour plots of vorticity 
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Figure 3.3: Mean and fundamental mode energies versus time for 346(0)5" wake. 
: 0.1 x mean flow energy. : Fundamental disturbance energy. 

at several times for the 346(0)5" run (an Re b = 346 wake with a two-dimensional 
fundamental only). Note that these plots are in a frame of reference moving with the 
free stream (moving left to right at speed U 0 > %,), hence the structures will appear 
to be convecting from right to left. The evolution starts with the linear growth 
of the disturbances, dominated by the growth of the most unstable (fundamental) 
wavelength perturbation (figure 3.2a). This is followed by nonlinear growth and 
rollup of the fundamental disturbance into a Karman vortex street (figure 3.2b), which 
consists of a staggered double row of roughly circular regions of vorticity, known as 
“rollers”. These are positive signed on the bottom of the wake, negative signed on the 
top, and separated by regions of low enstrophy density (figure 3.2e). In the absence of 
longer wavelength disturbances, the rollers become roughly circular in cross section 
(figure 3.2f) and settle into a stable configuration. The wake then slowly diffuses 
under the influence of viscosity. 

The process of linear growth, saturation, and decay described above is readily 
evident in figure 3.3 which shows the time evolution of the mode energy (the en- 
ergy contained in motions with a given streamwise wavelength) of the mean flow 
and fundamental disturbance Note that, in the figures, the mean energy line has 
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been rescaled to fit on the same plot as the fundamental mode energy. The funda- 
mental mode energy grows exponentially (linear growth region) until approximately 
t = 30, the fundamental mode then enters a nonlinear growth regime, and saturates 
at approximately t = 50 (at which time the vortex street has formed). Up to this 
time the fundamental mode has been drawing energy from the mean flow. The fun- 
damental mode then begins to decay, releasing energy into both shorter wavelength 
disturbances and back into the mean flow. Oscillations in both the fundamental and 
mean mode energies past a time of 100 represent a long period exchange of energy 
between the mean flow and fundamental. These oscillations damp out as the wake 
slowly approaches an asymptotic state in which the wake is dominated by diffusion 
of the vortex street. These oscillations will appear again in plots of the evolution of 
the halfwidth of this wake presented later in this chapter (figures 3.17 and 3.19). 

Returning to the enstrophy density plots, figure 3.2c and figure 3. 2d show a pinch- 
ing off event in the development of the 346(0)^ w &ke. Fluid containing vorticity of 
the opposite sign is convected across the wake centerline to the opposite side of the 
wake during the rollup process. This sort of event is not unique to the temporal wake. 
A similar convection of fluid across the wake is evident in the flow visualizations of the 
wake behind a circular cylinder by Zdravkovitch [59] and direct numerical simulations 
of a spatially evolving two-dimensional wake by Maekawa & Mansour [24]. 

Note that there is a slight asymmetry across the wake at the latest time shown 
in this simulation. This asymmetry is due to the accumulation of small errors intro- 
duced by the approximate dealiasing scheme, which introduces a very small amplitude 
random forcing at the longest wavelengths. Such small numerical errors are sufficient 
to break the exact numerical symmetry over a very long calculation, but do not have 
a significant impact on the overall results at the times examined. 

3.3.2 Effect of Disturbance Wavelength and Phase 

The presence of a fundamental wavelength disturbance causes the initially uniform 
wake to develop into the familiar double row of vortices that make up the Karman 
vortex street. The rows of vortices are staggered, with the array of rollers on one 
side of the wake being 180° out of phase with the the rollers on the other side. The 






CHAPTER 3. TWO-DIMENSIONAL SIMULATIONS 


59 



Figure 3.4: Iso-vorticity contours for 346(0)^ wake (two-dimensional fundamental 
and subharmonic) at various times, (a) t = 27.5. (b) t = 54.5. (c) t = 86.5. (d) t = 
119.3. (e) t - 144.7. (f) t = 200.8. (g) t = 304.1. Contours are 0.01 < |u>| < 0.4 in 
increments of 0.05. 
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(g) 


Figure 3.5: Iso-vorticity contours for 346(0)£|* wake (two-dimensional fundamental 
and subharmonic at at various times, (a) t = 27.6. (b) t = 54.7. (c) t = 85.6. 
(d) t = 118.7. (e) t = 144.6. (f) t = 157.1. (g) t = 198.7. Contours are 0.01 < |u>| < 
0.4 in increments of 0.05. 
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Figure 3.6: Iso-vorticity contours for 346(0)5?* wake (two-dimensional fundamental 
and subharmonic at |) at various times, (a) t = 27.6. (b) t = 54.6. (c) t = 84.9. 
(d) t = 107.2. (e) t = 119.6. (f) t = 155.4. (g) t = 247.7. (h) t = 301.1. Contours 
are 0.01 < |w| < 0.4 in increments of 0.05. 
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Figure 3.7: Mean, fundamental, and subharmonic mode energies versus time for 

346(0)^ wake. : 0.1 x mean flow energy. : Fundamental mode energy. 

: Subharmonic mode energy. 



Figure 3.8: Mean, fundamental, and subharmonic mode energies versus time for 

346(0)oi^. wake. : 0.1 x mean flow energy. : Fundamental mode energy. 

: Subharmonic mode energy. 
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Figure 3.9: Mean, fundamental, and subharmonic mode energies versus time for 

346(0)££* wake. : 0.1 x mean flow energy. : Fundamental mode energy. 

: Subharmonic mode energy. 

addition of longer wavelength disturbances has the effect of strengthening some of 
the rollers while simultaneously weakening others. This causes the rollers to merge, 
orbit one another, or shred as the wake evolves. The details of the evolution depend 
on which rollers are strengthened or weakened, and hence depends on the relative 
phasing of the fundamental and subharmonic disturbances. Figures 3.4, 3.5, and 
3.6 show the evolution of the vorticity fields for the wakes 346(0 )qq*, 346(0)55.*, and 
346(0)oix- As stated above, the initial development of all of these wakes parallels that 
of the flow in figure 3.2, which has only the fundamental wavelength disturbance, up 
through the development of the vortex street. Significant differences arise only after 
the fundamental disturbance begins to saturate during rollup. 

In the wake in figure 3.4 (346(0)5£*), the subharmonic, which is at zero phase with 
respect to the fundamental, acts to strengthen one of the rollers on the top of the wake 
while weakening the other (figure 3.4b). Because of the phasing, the two rollers on 
the bottom of the wake are initially of equal strength (the subharmonic disturbance 
has zero amplitude at the locations where the bottom rollers form). The strong top 
roller captures both of the bottom rollers (figure 3.4c), stretching them and forcing 
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them together (figure 3.4d). The elongated bottom rollers then agglomerate leaving a 
single new bottom roller and a pair of top rollers which are of roughly equal strength 
(figure 3.4e,f). One of the top rollers is then pulled to the bottom of the wake and 
sandwiched between the (periodic) bottom rollers. There it is stretched and slowly 
absorbed by the bottom roller (figure 3.4g) 

The wake in figure 3.5 (346(0)5|*), which has a subharmonic at a phase of tt/ 4 
with respect to the fundamental (shifted by one eighth of a fundamental wavelength), 
follows a similar evolution as the wake in figure 3.5. A dominant top roller forms which 
causes the bottom rollers, which in this case are of slightly different strengths, to 
agglomerate into a single large roller. This is followed by the bottom roller capturing 
one of the top rollers and forming a vortex pair. It is not clear whether this is a stable 
arrangement, but it persists well past t = 200. 

Simulations by Aref k Siggia [1] of a two-dimensional inviscid wake using a discrete 
vortex method showed similar results. When they initiated their simulations with a 
random variation in the position of their vortices (which translates into a random 
initial disturbance field) they observed pairing of regions of vorticity of the same sign 
and the formation of vortex pairs (as seen here) and triplets. The vortex pairs were 
observed to persist through the length of the simulations. 

Williamson k Roshko [56], in experiments on vortex formation in the wake of an 
oscillating cylinder, observed similar asymmetric pairing of vortices. By varying the 
frequency and amplitude of oscillation, they were able to produce repeatable wakes 
with various combinations of single vortices and pairs of opposite signed vortices. 
Under certain conditions the vortex pairs would self-convect away from the wake 
centerline (cf. their figure 17), much as is seen here in figures 3.4 and 3.5. 

The wake in figure 3.6 (346(0)ol^.), which has a subharmonic at a phase of tt/ 2 
with respect to the fundamental (shifted by one quarter of a fundamental wavelength), 
undergoes a markedly different evolution. In this simulation, the phasing of the 
subharmonic disturbance creates a symmetry condition such that each side of the 
wake develops a strong roller and a weak roller, with the difference in strength between 
the rollers on each side of the wake the same. Instead of one side of the wake pairing 
and then absorbing a roller from the other side, the symmetry of the flow allows the 
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rollers on both sides of the wake to pair simultaneously. Note that this results in 
a much faster growth of the new' vortex street than either of the other two wakes 
presented ( t ~ 160.0 as opposed to t > 200.0 for both 346(0)5ox an< ^ 346(0)oi* ). The 
resulting wake is nearly symmetric across the wake centerline (the small asymmetry 
is due to the fact that the value of 7r could only be represented to a finite accuracy), 
which is not the case for the wakes in figures 3.4 and 3.5. 

Maekawa, Mansour, k Buell [24] performed direct numerical simulations of the 
two-dimensional spatially evolving incompressible wake forced with high amplitude 
random noise. Groups of structures very similar to those which appear in figures 
3.4, 3.5, and 3.6 appear in those simulations (cf. their figure 11). At least for some 
structures in the spatially evolving wake, the random noise behaves dynamically like 
a fundamental plus a subharmonic at a random phase. 

Figures 3.7, 3.8, and 3.9 show the evolution of the mean, fundamental, and sub- 
harmonic mode energies for the wakes shown in figures 3.4, 3.5, and 3.6. Up to a time 
of approximately t — 30 the energies of the mean and fundamental mode develop as 
if there were no subharmonic disturbance present. This is consistent with the fact 
that up to that time the wake is in a linear growth regime, and therefore the main 
flow of energy should be from the mean to the fundamental. 

After a time of t = 30, the subharmonic has gained enough energy to begin to 
affect the development of the fundamental mode. The presence of the subharmonic 
disturbance, which initially saps energy from the fundamental, lowers the saturation 
(peak) energy of the fundamental disturbance and causes it to saturate slightly earlier. 
The subharmonic experiences a “mini-saturation” slightly after the fundamental mode 
saturates, and begins to decay until a time of approximately t = 60. The subharmonic 
then begins to rapidly extract energy from the fundamental disturbance. 

Up to a time of t = 90 the energies in the flows with a subharmonic disturbance 
develop identically, independent of phase. The time t = 90 corresponds to the time 
at which the rollers in the initial fundamental wavelength vortex street begin to pair 
(figures 3.4c/d, 3.5c/d, and 3.6c/d). After t = 90 evolution of the mode energies be- 
gins to become phase dependent. Between t = 90 and t = 150 all of the subharmonic 
wakes show a periodic exchange of energy between the fundamental and subharmonic 
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Figure 3.11: Mean, fundamental, subharmonic, and sub-subharmonic mode energies 

versus time for 346(0 )^q wake. : 0.1 x mean flow energy. : Fundamental 

mode energy. : Subharmonic mode energy. : Sub-subharmonic mode 

energy. 

disturbances, as well as between the fundamental and the mean (this is most obvious 
in figure 3.9, but appears to a lesser degree in the other two figures). In each case, 
the fundamental eventually decays, dumping its energy into both the subharmonic 
and the mean, leaving the subharmonic as the dominant disturbance mode. 

All of the wakes with subharmonic disturbances develop large structures that 
could potentially be constricted by the finite streamwise size of the box. To attempt 
to get a handle on the significance this constriction on the development of the wakes 
in figures 3.4 through 3.6, a wake with a sub-subharmonic disturbance was simulated. 
Figure 3.10 shows plots of iso-vorticity contours at four times in the development of 
the 346(0)ggo wake. Comparing figure 3.10 with the results for the 346(0)gJ* wake in 
figure 3.4, it is clear that at least up to a time of t = 200, the finite length of the 
computational box is not unduly affecting the wake structure. The sub-subharmonic 
disturbance has only a small noticeable effect on the development of the wake. Only at 
very late times is the presence of the sub-subharmonic disturbance likely to become 
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important, but for the moderate times studied here and in the three-dimensional 
simulations it is reasonable to assume that a subharmonic length box is sufficient. 

Figure 3.11 shows the evolution of the mean, fundamental, subharmonic, and 
sub-subharmonic mode energies for the wake in figure 3.10. Comparing this figure 
to the corresponding figure for the 346(0 )qoi wake (figure 3.7 on page 66) it is clear 
that the sub-subharmonic disturbance remains at a very low energy late until in the 
simulation. The sub-subharmonic has no noticeable impact on the evolution of the 
energies of the shorter wavelength disturbances up until a time of approximate!} 
t = 150. After that time, the sub-subharmonic disturbance begins to slowly draw 
energy from the (dominant) subharmonic mode. Even at the latest time simulated 
the sub-subharmonic is beginning to grow rapidly at the expense of the subharmonic, 
but is still much weaker than the subharmonic mode. 

3.3.3 Effect of Reynolds Number 

Figures 3.12, 3.13, 3.14, and 3.15 show contour plots of vorticity for sets of wakes 
with the same initial disturbances at approximately the same time, but with different 
Reynolds numbers. In general, varying the Reynolds number has minimal effect on the 
large scale structure of the flow (given a high enough Reynolds number to allow the 
initial rollup). The major differences between flows with different Reynolds numbers 
has to do with the development of small scales. 

At very low Reynolds numbers (figure 3.12a), viscosity dominates the wake to the 
point that it never enters the nonlinear growth regime. At slightly higher Reynolds 
numbers (figures 3.12b, 3.13a 3.14a, and 3.15a), the wake rolls up and non-linear 
dynamics, including pairing of the rollers, are evident. The structures are still very 
diffuse, however, and adjacent regions of opposite signed vorticity annihilate each 
other quickly. 

As the Reynolds number increases (figures 3.12c-f, 3.13b,c 3.14b,c, and 3.15b,c), 
the decay of the intensity of the large vortical structures is significantly reduced 
and much smaller vorticity scales begin to appear. These small scales are the result 
of vortical regions being wrapped and folded around the rollers and vorticity being 
stretched in the high strain-rate regions between rollers. Without the strong diffusion 
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Figure 3.12: Iso-vorticity contours for (?)(0)J” wakes (two-dimensional fundamental 
with various Reynolds numbers), (a) Reb = 69, t = 200.7. (b) Reb = 119 , t = 201.8. 
(c) Re b = 346, t = 198.6. (d) Re b = 692, t = 201.3. (e) Re b = 1384, t = 197.8. 
(f) Re b = 2768, t = 196.3. Contours are 0.01 < |w| < 0.4 in increments of 0.05. 
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(c) (d) 


Figure 3.14: Iso-vorticity contours for (?)(0)£|* wakes (two-dimensional fundamental 
plus subharmonic at j with various Reynolds numbers), (a) Ret, = 119, t = 201.4. 
(b) Re b = 346, t = 198.7. (c) Re b = 692, t = 193.4. (d) Re b = 1384, t = 195.2. 
Contours are 0.01 < |u?| < 0.4 in increments of 0.05. Arrow indicates rollup of 
secondary shear layer. 
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Figure 3.15: Iso-vorticity contours for (?)(0)oi* wakes (two-dimensional fundamental 
plus subharmonic at | with various Reynolds numbers), (a) Re b = 119, t = 205.3. 
(b) Re b = 346, t = 196.0. (c) Re b = 692, t = 197.3. (d) Re h = 1384, t = 200.2. 
Contours are 0.01 < |u>| < 0.4 in increments of 0.05. 
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present in the wakes with lower Reynolds number, these fine structures do not get 
smeared out or absorbed by the larger structures. 

As the Reynolds number increases there is also an increase in the strength and 
persistence of the vorticity in the fluid which has been convected across the wake 
in pinching off events during rollup and pairing. These pockets of pinched off fluid 
contain vorticity of the opposite sign as the rollers, adding significant complexity to 
the overall wake structure. 

For the cases with pairing, at the higher Reynolds numbers (figures 3.13c,d 3.14c,d, 
and 3.15c,d) the post-pairing rollers show small scale internal structures (this is most 
evident in figure 3.13d, where the bottom “roller” still consists of two separate struc- 
tures). There is also evidence of the onset of secondary rollup of the small scale shear 
layers which develop during the pairing process. The most obvious of these is in the 
long thin structure at the top of figure 3.14d (see arrow). 

3.3.4 Scalar Field 

Figure 3.16 shows contours of enstrophy density and scalar concentration for three of 
the wakes discussed above. Since the initial scalar field w r as set to the initial enstrophy 
density field, it is reasonable to assume that contour plots of these quantities should 
be similar in appearance. Any major differences are due to the fact that vorticity 
can be canceled out by vorticity of the opposite sign whereas the scalar is a positive 
quantity which is conserved. 

It is clear from figure 3.16 that even at high Reynolds numbers where there is 
complex structure to the wake, the scalar field and the enstrophy density field are 
similar. Portions of the flow which have large enstrophy density tend to be regions 
where vorticity of opposite signs is not present. These regions have matching large 
values for the passive scalar. Regions of the flow with nonzero passive scalar con- 
centrations but little enstrophy density correspond to once-vortical fluid which has 
undergone cancellation. 

The scalar concentration plots in figure 3.16 give a good estimate for the level 
of entrainment of free-stream fluid by the various wakes. The wake in figure 3.16b 
(1384(0)oo^) has entrained a large quantity of fluid from both the freestream above 
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(e) (f) 


Figure 3.16: Comparison of iso-enstrophy density and iso-scalar contours for 

1384(0)q^ x wakes (two-dimensional fundamental and subharmonic at various phases). 
1384(0 )oq£\ t = 208.1. (a) Enstrophy Density, (b) Passive scalar. 1384(0 )q£ x . 

t = 195.2. (c) Enstrophy Density, (d) Passive scalar. 1384(0 )q£ x , t = 200.2. (e) En- 
strophy Density, (f) Passive scalar. Contours are 0.05 < M,Ci < 0.4 in increments 
of 0.05. Arrows indicate entrained fluid. 
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and below the wake (see arrows). In contrast, the wake in figure 3. 1 6f (1384(0 )q|^.) has 
entrained very little freestream fluid (only what has been convected in between the 
periodic rollers on each side of the wake — see arrows). The differences in entrainment 
make a significant difference in the mixing which occurs in these wakes, and has a 
strong impact on their relative growth rates. 

3.4 Growth of the Mean Flow 

In all of the simulations, the mean velocity profile is defined to be the average velocity 
over the streamwise and spanwise directions for a given cross-stream position. The 
(“direct”) wake halfwidth, 6, is then calculated by finding the maximum mean velocity 
for a given time in the simulation and taking half the width between outermost 
crossings of the mean velocity profile with the half maximum velocity point. 

Since many of the flows to be presented develop non-Gaussian mean profiles, 
primarily due to the limited number of long wavelength modes available over which 
to average, a second, integral measure will also be presented. This integral halfwidth, 
6,, is defined as 

b,= 

w’here A u is the mean (streamw’ise) velocity defect profile. The integral halfwidth 6, 
has been formulated such that for a Gaussian velocity defect profile 6, = b. Since the 
initial velocity profiles are all Gaussian, b t (t = 0) = bo = 1.0. 

3.4.1 Effect of Disturbance Wavelength 

Similarity arguments applied to the temporally developing plane wake, where the 
momentum per unit streamwbse length is assumed to be the only important governing 
parameter, lead to the conclusion that the far wake should grow asymptotically like 
t*. Therefore it is convenient to use plots of the square of the wake halfwidth, b 2 or 
b ■? versus time, t, to study the growth of the mean flow. In these coordinates, b ~ t * 
or b x ~ t 2 growth will appear as a straight line. 
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(a) 



(b) 

Figure 3.17: Square of normalized wake halfwidth versus normalized time for Re = 
346 wakes with various disturbance wavelengths, o: 346(0) O: 346(0)££. □: 
346(0 )qJo • ( a ) Direct halfwidth, b. (b) Integral halfwidth b % . 
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As a reference, the streamwise fundamental wavelength, as calculated from linear 
stability theory, is Lj = 7.85. The halfwidth b is therefore equal to the halflength of 
the computational domain when b 2 = 15.4 for the fundamental only case, b 2 = 61.7 for 
the cases with subharmonic disturbances, and b 2 = 246.7 for the cases with sub-sub- 
harmonic disturbances. When b 2 approaches these values it is likely that the stream- 
wise confinement imposed by the computational box will be become verj< significant. 
(Note: these values will also hold for the later three-dimensional computations). 

Figure 3.17 shows plots of the square of the wake halfwidths versus time for 
three wakes with i?ej, = 346 started from initial fields containing a fundamental, 
a fundamental and a subharmonic, and a fundamental, a subharmonic, and a sub 
subharmonic respectively. The circles correspond to the simulation that appears in 
figure 3.2 (346(0)g“), the diamonds to the simulation in figure 3.4 (346(0)g^), and 
the squares to the simulation in figure 3.10 (346(0);^). 

As discussed above, the wakes go through three distinct stages of evolution. Up 
to a time of approximately t = 30, which is the linear growth regime, the flow is 
laminar with exponentially growing sinusoidal disturbances. After a time of t 30, 
the flow enters a nonlinear growth regime where the fundamental disturbance grows 
rapidly and the wake spreads quickly. At a time of approximately t = 50 the wake 
width peaks as the fundamental begins to saturate and the vorticity field rolls up to 
form the Karman vortex street. As the vortex street forms, the rollers relax slightly 
toward the wake centerline and the wake halfwidth drops. 

Up to a time of t = 90, the growth of the wakes shown (and indeed all of the tw’o- 
dimensional wakes studied) is clearly dominated by the evolution of the fundamental 
disturbance. After a time of t = 100, however, the longer wavelength disturbances 
become important. After t = 100, the simulation with only the fundamental distur- 
bance (346(0)5**) settles down to a uniform vortex street which spreads only slowly 
through viscous diffusion. The small oscillations in wake width at late time are due to 
the vortices positioning themselves into an asymptotically stable configuration, and 
are therefore less pronounced in the plot of the integral halfwidth. 
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The sharp increase in wake width at a time of / = 105 for the simulations with the 
longer wavelength (subharmonic, 346(0)^, and sub-subharmonic, 346(0)5 qq) distur- 
bances corresponds to the growth of the subharmonic disturbance, which results in 
the initial pairing of the vortices on the bottom side of the wake (figures 3.4c,d). The 
wake width then peaks as the bottom rollers finish pairing and force one of the top 
rollers upwards (figures 3.4e). As the wake approaches a subharmonic vortex street 
configuration the rollers again settle back toward the centerline and the wake width 
decreases. 

There is very little difference in the evolution of the wake with only the subhar- 
monic and the wake with the addition of a sub-subharmonic disturbance. Up to the 
time simulated, the sub-subharmonic does not play a significant role. For a very long 
simulation, however, it is reasonable to expect that the sub-subharmonic would even- 
tually lead to a second pairing of the vortex street. The effect of the sub-subharmonic 
are more evident in the plots of the integral halfwidth at late times, but the difference 
is still not significant. 

3.4.2 Effect of Disturbance Phase 

Figure 3.18 shows plots of the wake halfwidths versus time for the three simulations 
with Reb = 346 which were initiated with a fundamental and a subharmonic distur- 
bance. These are the same wakes that appear in the vorticity contour plots in figures 
3.4, 3.5, and 3.6. As one would expect after examining the vorticity contour plots, the 
phasing of the subharmonic has a very strong impact on the evolution of the mean 
width. 

The two wakes with subharmonic phasing such that the rollers pair on one side and 
then capture a roller from the other side (346(0)5oi and 346(0)5§.* ) have significantly 
greater maximum widths than the wake that pairs on both sides simultaneously. The 
point of maximum width for all of the wakes occurs in the middle of the pairing 
process, before the original rollers have fully combined. In the wakes that pair only 
on one side, the pairing process pushes one of the rollers on the other side of the wake 
away from the centerline. This serves to sharply increase the mean halfwidth of the 
wake. After the pairing rollers have completed their agglomeration, the roller that 
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( a ) 



(b) 

Figure 3.18: Square of normalized wake halfwidth versus normalized time for 

346(0)jvfw wakes (two-dimensional fundamental and subharmonic at various phases), 
o: 346(0)go*. O: 346(0)g|*. □: 346(0)gi*. (a) Direct halfwidth 6. (b) Integral 
halfwidth b % . 
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was pushed out is drawn back toward the centerline, decreasing the halfwidth again. 
In the wake which pairs on both sides simultaneously, the rollers on both sides of the 
wake stay close to the centerline, and the width peaks only slightly. 

3.4.3 Effect of Reynolds Number 

Figures 3.19, 3.20, 3.21, and 3.22 show plots of the square of the halfwidths versus time 
for wakes with identical initial disturbances but with different Reynolds numbers. The 
wake with the lowest Reynolds number shown is 69(0)5** (°P en diamonds in figure 
3.19). At this Reynolds number, the wake never rolls up. Although there is some 
nonlinear growth of the disturbances at early times (as evidenced by the small hump 
in the growth curves), the growth of this wake is due almost entirely to (rather rapid) 
viscous diffusion. 

At a somewhat higher Reynolds number (119(0)5**, 119(0)55*, 119(0)5|*, and 
119(0)5|*, all shown as open squares in figures 3.19, 3.20, 3.21, and 3.22), the wakes 
roll up to form vortex streets. Although pairing does occur where a subharmonic 
disturbance is present, the late time growth of all the wakes quickly become dominated 
by diffusion. 

As the Reynolds number increases further, some measure of Reynolds number 
independence begins to appear in the growth of the mean flow. This is particularly 
true for the (?)(0)5** wakes (figure 3.19). The growth curves for the Rtb = 346. 
692, 1384, and 2768 wakes lie almost exactly atop one another. The only noticeable 
difference is that the initial wake rollup occurs slightly earlier at higher Reynolds 
numbers. 

For the wakes that undergo pairing (figures 3.20, 3.21, and 3.22) the case for 
Reynolds number independence is somewhat weaker. The growth seems to be inde- 
pendent of Reynolds number up to approximately midway through the pairing of the 
rollers, at a point where the rollers are still recognizable as separate structures. Be- 
yond this point the growth of the halfwidth appears to become dependent on viscosity 
again. 

A large part of the variation in the “direct” halfwidth 6, as compared to the integral 
halfwidth 6,, is due to the sensitivity of the measure to the shape of the mean velocity 
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(a) 



(b) 

Figure 3.19: Square of normalized wake halfwidth versus normalized time for (?)(0)5" 
wakes (two-dimensional fundamental with various Reynolds numbers). O: 69(0)5". 
□: 11 9(0)5“. o: 346(0)5". A: 692(0)5". V= 1384(0)5". x: 2768(0)5". (a) Direct 
halfwidth b. (b) Integral halfwidth 
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(a) 



(b) 

Figure 3.20: Square of normalized wake halfwidth versus normalized time for (?)(0)g^ 
wakes (two-dimensional fundamental and subharmonic with various Reynolds num- 
bers). □ : 119(0)S£. o: 346(0)5^. A: 692(0)g£. 1384(0)££. (a) Direct halfwidth 

b. (b) Integral halfwidth 6,. 
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(b) 

Figure 3.21 : Square of normalized wake halfwidth versus normalized time for (?)(0)5|* 
wakes (two-dimensional fundamental and subharmonic at ^ with various Reynolds 
numbers). □: 119(0)$f*. o: 346(0)g|*. A: 692(0)S|*. V= 1384(0)gi*. (a) Direct 

halfwidth b. (b) Integral halfwidth b t . 
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(b) 

Figure 3.22: Square of normalized wake halfwidth versus normalized time for (?)(0)g£* 
wakes (two-dimensional fundamental and subharmonic at | with various Reynolds 
numbers). □: 119(0)g| r r . o: 346(0)g|*. A: 692(0)51*. V : 1384(0)£|*. (a) Direct 
halfwidth 6. (b) Integral halfwidth b t . 
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Figure 3.23: Comparison of mean velocity profiles for (692)(0)5o£ and (1384}(0)gi£ 
wakes (two-dimensional fundamental and subharmonic at 0). A: 692(0)qox > ^ = 199.1. 
y: 1384(0 )qox ? t = 208.1. (Symbols correspond to those in figure 3.20). Arrows 
indicate locations for measurement of “direct” halfwidths. 

profile. The “direct” halfwidth is simply taken as half the width between the two 
outermost points where the mean streamwise velocity is at the half maximum velocity 
defect level. If the mean defect velocity profile is non-Gaussian, and in particular if 
it has “shoulders” or multiple peaks, then large variations in the measured halfwidth 
can appear between flows -.bat look very similar in terms of vorticity distribution. 
This is an inherent shortcoming of the “direct” halfwidth measure b which becomes 
particularly acute in cases such as this where the mean profiles can be far from 
Gaussian. 

This is illustrated by figure 3.23, which shows mean streamwise velocity profiles 
for the simulations that are presented in figures 3.13c and 3.13d. The vorticity fields 
of these flows appear very similar, but due to the way the large scale structures 
have arranged themselves, they produce very different mean profiles. At the lower 
Reynolds number, the rollers are more scattered and diffuse, which results in a broad, 
flattened mean velocity profile. At the higher Reynolds number, however, the fact 
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that the rollers are compact, and that the top roller has been pulled further toward 
the bottom of the wake produces a mean profile with a narrow', intense peak. Though 
the vorticity in both wakes has spread out over approximately the same spanwise 
extent (which is another possible measure of the width), the “direct” halfwidth b 
(arrows) is very different between the two cases. 

As noted above, the difference in the integral halfwidth, 6^, is much smaller than 
for the “direct” halfwidth b. There is, however, still a significant difference in the 
widths of the 692(0 )qo£ and 1384(0)gj£ wakes (A and V in figure 3.20b respectively). 
The higher Reynolds number 1384(0)^ wake is significantly more compact than 
the lower Reynolds number flows. This is due to the fact that at this particular 
subharmonic phasing, one side the wake develops a single region of concentrated 
vorticity, while the vorticity in the other side of the w'ake remains in two weaker 
regions. The difference in the relative strength of the regions of vorticity that form 
allow the single strong region to capture one of both of the weaker regions from the 
opposite side of the wake. At the highest Reynolds number shown, 1384(0)^, the 
very strong concentrated vorticity region at the bottom of the wake captures both 
of the weaker vorticity concentrations at the top of the wake, resulting in a very 
compact configuration. At the lower Reynolds numbers shown just one of the upper 
vorticity concentrations is captured, resulting in a wider wake (see figure 3.13 on 
page 75). In contrast, the other two subharmonic phasings examined result in wakes 
in which the vorticity concentrations that develop on opposite sides of the wake are 
more balanced in strength. Neither side of the wake can capture all of the significant 
vorticity concentrations from the other side of the wake, independent of the wake 
Reynolds number (see figures 3.14 and 3.15 on pages 76 and 77). 



Chapter 4 

Three-dimensional Simulations 


4.1 Motivation 

To explore the effects of initial conditions on the development of the three-dimen- 
sional incompressible plane wake, a set of three-dimensional simulations has been run 
using various initial conditions. The effects of disturbance wavelength, phasing, and 
wake Reynolds number are examined in this section. 


4.2 Simulation Parameters 

4.2.1 Three-dimensional Forcing 

As discussed in section 2.7, the initial conditions for the three-dimensional simulations 
were composed of a Gaussian mean streamwise velocity profile with sets of small am- 
plitude periodic disturbances superimposed. Simulations were performed at Reynolds 
numbers based on initial halfwidth and defect velocity of between 69 and 2768. The 
majority of wakes simulated were at a Reynolds number of Re^ — 346. Based on the 
results of the two-dimensional simulations presented in chapter 3, this should be a 
sufficiently high Reynolds number to capture the main effects of the variety of initial 
disturbances on the large scale development of the wake. Simulations at the highest 
Reynolds numbers, which are much more computationally demanding, were run for 
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only a few select sets of initial disturbances in order to examine the transition to 
turbulence. Table H.2 in appendix H gives a summary of all of the three-dimensional 
simulations that have been performed. 


4.3 Structural Development of the Three-dimen- 
sional Wake 

Figure 4.1 shows perspective views of three-dimensional iso-enstrophy density con- 
tours at four times in the development of the 346(60)g°* wake (two-dimensional fun- 
damental disturbance and a pair of 60° oblique disturbances at the subharmonic 
wavelength). The freestream flow direction is from left to right. The contour level 
was chosen to be |w| = 0.4. This corresponds to approximately 60% of the initial 
mean wake defect velocity divided by the initial wake halfwidth, which serves as a 
rough measure of the mean gradient in the initial flow. 

The development of this wake is typical of the wakes initiated with a two-dim- 
ensional fundamental disturbance and a pair of 60° oblique disturbances at the sub- 
harmonic wavelength. The initial development is primarily two-dimensional, with 
the fundamental (Karman) mode growing most quickly as predicted by linear theory. 
The flow develops well defined spanwise rollers with relatively weak streamwise struc- 
tures stretching between rollers on the same side of the wake (figure 4.1b). Once the 
rollers are established, the three dimensional disturbances grow rapidly in strength 
under the influence of the resulting straining field, (figure 4.1b,c). As the stream- 
wise structures become more intense, they begin to distort the rollers, breaking up 
their spanwise coherence (figure 4.1c). Eventually, the streamwise structures, both 
the original structures and ones which are the result of the distortion of the spanwise 
rollers, become the dominant features in the flow (figure 4. Id). Note that as the wake 
spreads with time, the maximum vorticity due to the mean flow drops, so the mag- 
nitude of a fixed level of vorticity relative to the gradients of the mean flow becomes 
substantially larger. The structures which appear in figure 4. Id are at an intensity 
level far above the vorticity due to the mean shear. Compared to the contour level of 
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(c) 


Figure 4.1: See caption page 97. 
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Figure 4.1: Iso-enstrophy density contour for 346(60)£°* wake (two-dimensional 

fundamental plus three-dimensional subharmonic). Contour level is |u.’| = 0.4. 
(a) t = 22.8. (b) t = 52.8. (c) t = 102.7. (d) t = 204.8. 
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Figure 4.2: Iso-enstrophy density contour for 346(60)g** wake (two-dimensional fun- 
damental plus three-dimensional fundamental), t = 196.9. Contour level is |u>| = 0.2. 

|o>| = 0.4 shown, the mean shear at that late time, based on halfwidth and maximum 
defect velocity, is 0.037, more than an order of magnitude lower. 

4.3.1 Effect of Disturbance Wavelength 

In stark contrast to the wake shown in figure 4.1 is the wake shown in figure 4.2. This 
is the 346(60)°** wake (which has a two-dimensional fundamental disturbance and a 
pair of oblique disturbances at the fundamental wavelength). The contour level here 
is 0.2, half the level used in figure 4.1. With these initial conditions, strong three- 
dimensionality fails to develop even by the advanced time shown. The streamwise 
structures are present (the dominant ones can be seen running between rollers on 
the bottom side of the wake at this disturbance phasing), but exist as broad flat 
regions of vorticity, canted to the plane of the wake, as opposed to the intense tube- 
like structures which appear in the 346(60)^* wake. The relatively low circulation 
of these weak streamwise structures is evidenced by their minimal impact on the 
spanwise rollers. The rollers in the wake in figure 4.2 are only slightly corrugated. 
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whereas at the corresponding time for the 346(60 )q°* wake (figure 4. Id), the spanwise 
rollers have been entirely overwhelmed by the streamwise structures. 

The effects of the addition of longer wavelength disturbances are illustrated by the 
wakes in figures 4.3, 4.5, and 4.7. Figure 4.3 shows the results from the 346(60 )qq£ sim- 
ulation. This wake was started from initial conditions similar to those used in figure 
4.2, but with the addition of a two-dimensional subharmonic disturbance. Again, the 
evolution is dominated by two-dimensional dynamics, even at the late time (t = 202.8) 
shown. As in the corresponding two-dimensional wake (figure 3.4), the two bottom 
rollers have paired to form one coherent roller, captured one of the top rollers, and 
drawn it down to the bottom of the wake. Although the three-dimensional structures 
are stronger than those found in the 346(60)^ simulation, and have a tube-like shape 
as opposed to a flat shape, they are still substantially weaker than those found in the 
346(60)g°* simulation. The effect of the moderate strength streamwise structures is 
to break up the spanwise coherence of the captured roller, with little impact on the 
other two (primary) rollers. 

The essential two-dimensionality of this flow is reinforced by the plot in figure 
4.4. This figure shows a comparison of the evolution of the first few streamwise 
(two-dimensional) mode energies for the three-dimensional wake 346(60)^ (symbols) 
and the corresponding two-dimensional wake 346(60)gjf (lines). Up to a time of 
approximately t = 150 there is no significant difference in the development of the 
mean, two-dimensional fundamental or subharmonic mode energies between the two- 
dimensional and three-dimensional wake. Even at very late times the differences are 
relatively small. The subharmonic mode energy for the three-dimensional wake is 
somewhat lower due to energy being transferred into the three-dimensional modes, 
but the wake dynamics are still essentially those of a two-dimensional wake. 

Figure 4.5 shows an iso-enstrophy density contour for the 346(60)^ simulation, 
where a two-dimensional subharmonic has been added to the initial disturbance field. 
Note that the contour level shown is again |u>| = 0.4. The addition of a two-dimen- 
sional subharmonic disturbance does not have a major impact on the development 
of the streamwise structures as compared to the 346(60)5°* simulation in figure 4.1. 
Similar primary streamwise structures exist in both flows, running between rollers 
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Figure 4.3: Iso-enstrophy density contour for 346(60)$£ wake (two-dimensional fun- 
damental and subharmonic plus three-dimensional fundamental), t = 202.8. Contour 
level is |u;| = 0.2. 
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Figure 4.4: Comparison of mean, two-dimensional fundamental, and two-dimensional 
subharmonic mode energies versus time for 346(0)55* (lines) and 346(0 )q 5* (symbols) 

wakes. , o: 0.1 x mean flow energy. , O: Fundamental mode energy. 

, □: Subharmonic mode energy. 


on the same side of the wake. The main effect of the two-dimensional subharmonic 
disturbance is on the spacing of the rollers, causing the two bottom rollers and one of 
the top rollers to cluster into a tight group and begin to pair as was observed in the 
corresponding two-dimensional wake in figure 3.4d and the weakly three-dimensional 
wake in figure 4.3. In this wake, however, the strong streamwise structures inhibit 
the pairing process by distorting the rollers. 

This is more easily seen by examining the evolution of the streamwise mode ener- 
gies. Figure 4.6 shows a comparison of the disturbance energies for the two-dimension- 
al 346(60)55£ wake (lines) and the three-dimensional 346(60)^ wake (symbols). Note 
that both the fundamental and subharmonic disturbances in the three-dimensional 
wake peak earlier and at a lower energy than in the corresponding two-dimensional 
wake. Furthermore, the late time energy of the subharmonic disturbance is signifi- 
cantly lower in the three-dimensional wake. This is a result of the three-dimensional 
motions drawing energy out of the long wavelength two-dimensional modes. 
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Figure 4.5: Iso-enstrophy density contour for 346(60)^ wake (two-dimensional fun- 
damental and subharmonic plus three-dimensional subharmonic), t = 96.7. Contour 
level is |u>| = 0.4. 
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Figure 4.6: Comparison of mean, two-dimensional fundamental, and two-dimensional 
subharmonic mode energies versus time for 346(0)^ (lines) and 346(0 )§q£ (symbols) 

wakes. , o: 0.1 x mean flow energy. , O: Fundamental mode energy. , 

□ : Subharmonic mode energy. 



Figure 4.7: Iso-enstrophy density contour for 346(60)5^ wake (two-dimensional fun- 
damental, subharmonic, and sub-subharmonic plus three-dimensional subharmonic 
and sub-subharmonic), t = 103.2. Contour level is |u>| = 0.4. 
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Figure 4.8: Mean and two-dimensional fundamental, two-dimensional subharmonic, 
and two-dimensional sub-subharmonic mode energies versus time for 346(0)5°° wake. 
: 0.1 x mean flow energy. : Fundamental mode energy. : Subhar- 
monic mode energy. : Sub-subharmonic mode energy. 

Figure 4.7 shows iso-enstrophy density contours for the 346(60)5oo wake, which 
has the longest wavelength disturbances studied. Note that the streamwise extent of 
the field shown in figure 4.7 is twice that of the fields in the previous figures (four 
fundamental wavelengths long instead of two). The effect of the additional two-dim- 
ensional sub-subharmonic disturbance is rather minimal up to the time shown. It has 
caused the sets of paired rollers to be offset slightly in the cross-stream direction, but 
there is no indication that a second pairing is likely. This conclusion is supported 
by figure 4.8 which shows the evolution of the energies of the first few streamwise 
disturbance modes. Note that the sub-subharmonic disturbance (dotted line) is the 
weakest of the disturbance modes shows, and is decaying rapidly at late times. Since 
the sub-subharmonic mode is required for a second pairing to occur, such a paring is 
impossible in the time frame observed. 

At late times the the spanwise structures are very similar to the structures in the 
wake in figure 4.5. The effect the addition of the oblique sub-subharmonic disturbance 
is to impart a spanwise variation in the strengths of both the streamwise structures 
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r — x 

Figure 4.9: Iso-enstrophy density contour for 346(60) Ox x wake (two-dimensional fun- 
damental plus three-dimensional subharmonic at p. t = 96.5. Contour level is 
M - 0.4. 

and the spanwise rollers. This small increase in the complexity of the overall flow has 
little impact on the development of the wake. 

4.3.2 Effect of Disturbance Phase 

Figures 4.9, 4.10, 4.11, and 4.12 illustrate the effects of phasing of the oblique distur- 
bance on the structure of the three-dimensional wakes discussed above. The general 
effect is to decrease the strength of the streamwise structures on one side of the wake 
and increase the strength of the streamwise structures on the other. For the wakes 
initiated with an oblique disturbance at the fundamental wavelength (figures 4.10 and 
4.2, and figures 4.11 and 4.3) this results in the already weak streamwise structures 
on the bottom of the wake dropping in intensity below the contour level shown in 
favor of streamwise structures on the top of the wake, which remain too weak to 
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Figure 4.10: Iso-enstrophy density contour for 346(60)0^ wake (two-dimensional 
fundamental plus three-dimensional fundamental at ^). t = 197.0. Contour level is 

M = 0.2. 

appear. For the wakes initiated with a oblique disturbance at the subharmonic wave- 
length (figures 4.9 and 4.1, and figures 4.12 and 4.5) this results in a balancing of the 
strengths of the streamwise structures, which previously favored the top of the wake. 

4.3.3 Effect of Reynolds Number 

Figure 4.13 shows a series of iso-enstrophy density contours for (?)(60)j^ wakes all 
at approximately the same time, but with progressively increasing Reynolds number. 
Note that as the Reynolds number increases, the contour level shown also increases. 

Even at a rather low Reynolds number (the 119(60)^ wake in figure 4.13a), 
the streamwise structures which arise from the oblique subharmonic disturbance are 
quite strong and dynamically important. Though they do not overwhelm the spanwise 
rollers, they do distort the rollers significantly. 

As the Reynolds number increases, the relative strengths of the streamwise struc- 
tures also increase, the spanwise rollers become more distorted, and finer flow scales 
appear. Even in the highest Reynolds number wake simulated (the 2768(60)5°* wake 
in figure 4.13e), the influence of the oblique subharmonic is apparent. Very intense 
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Figure 4.11: Iso-enstrophy density contour for 346(60) ( ^ 1 wake (two-dimensional 
fundamental and subharmonic plus three-dimensional fundamental at -j). t = 192.6. 
Contour level is |u>| = 0.2. 
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Figure 4.12: Iso-enstrophv density contour for 346(60) O o 2 x wake (two-dimensional 
fundamental and subharmonic plus three-dimensional subharmonic at |). t = 97.3. 
Contour level is |u?| = 0.4. 
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(a) 



(b) 


Figure 4.13: See caption page 111. 
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Figure 4.13: Iso-enstrophy density contours for (?)(60)j££ wakes (two-dimension- 
al fundamental plus three-dimensional subharmonic) at various Reynolds numbers, 
(a) 119(60)5°*, t = 97.9. \u>\ = 0.2. (b) 346(60)g£, t = 102.7. M = 0.4. 
(c) 692(60)§°5, t = 101.2. |w| = 0.8. (d) 1384(60)g£, t = 101.5. M = 1.6. 
(e) 2768(60)5°^, t = 102.4. |w| = 2.4. 
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Figure 4.14: Iso-enstrophy density contour for 1384(60)°** wake (two-dimensional 
fundamental plus three-dimensional fundamental), t = 100.2. Contour level is |u?| = 
0.2. 

fine scale motions occur at locations where the coherent streamwise structures appear 
in the lower Reynolds number flows. This is an indication that the generation of small 
scales in the wake is linked the very intense straining fields created by the streamwise 
structures. 

Though the early development of all the wakes calculated is similar, the devel- 
opment of significant strong three-dimensionality does not occur in the absence of 
the oblique subharmonic. Figure 4.14 shows an iso-enstrophy density contour for 
the 1384(60)oU wake. The streamwise structures are evident, and, as evidenced by 
the visible corrugation in the spanwise rollers, have a somewhat greater impact on 
the flow dynamics than in the corresponding lower Reynolds number (?)(60)[J** wake 
presented previously. However, the three-dimensional motions are still very weak, 
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and the wake remains nearly two dimensional. Similar results are found for computa- 
tions started from the same set of initial disturbances, but with different disturbance 
phasings. 

This result differs markedly from the compressible wake computed by Chen et 
ai [10], which showed strong three-dimensional development for certain disturbance 
phases and only weak three-dimensional development for others. It is unclear why 
this difference exists. Chen et ai. attributed the sensitivity to three-dimension- 
al disturbance phase in the compressible wake to the fact that at certain phasings 
significant amounts of vorticity from the three-dimensional mode resides in the high 
strain-rate regions between the spanwise rollers where it can be amplified through 
stretching. Other phases are such that most of the vorticity from the three-dimen- 
sional mode gets wrapped into the relative low strain-rate rollers, where it is only 
weakly amplified. If this were the only process occurring the same results should hold 
for the incompressible wake. A more careful comparison of the respective datasets 
needs to made before the reason for the observed difference can be determined. 

At this point it is prudent to check that the wake simulations presented, particu- 
larly the simulations of the higher Reynolds number wakes, are well resolved. To do 
this it should be sufficient to demonstrate that the most extreme case, the Re b = 2768 
wake, is adequately resolved. This was accomplished by calculating the Taylor mi- 
croscales as was done for the test case in section 2.8.3 (see appendix A for definitions 
of the microscales used). Analysis of the 2768(60)*!* wake yielded a minimum (as a 
function of x 2 ) Taylor microscales of Aj = 0.26 and Xti = 0.35 respectively. The least 
well resolved direction in this flow is the cross-stream, with a grid size of Ax 2 = 0.085 
which translates to over three grid cells per microscale length. This resolution is 
adequate. 

This conclusion is reinforced by by examining the energy spectra of the simulated 
wakes. Figure 4.15a shows the energy spectra for the 346(60)jg* wake and figure 
4.15b shows the energy spectra for the 2768(60);^ wake. The 2768(60)** wake is 
the most demanding of the simulations performed. Since the spectra of both wakes 
show large, clean dropoff of their spectra at high wavenumbers it is clear that both 
simulations are sufficiently resolved. 
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(a) 



(b) 


Figure 4.15: Streamwise energy spectra for 346(60)^ and 2768(60)5°* wakes. : 

E(k i). (a) 346(60)^, t = 102.7. (b) 2768(60)*!*, t = 102.4. Arrows indicate 
minimum grid wavenumber. 





CHAPTER 4. THREE-DIMENSIONAL SIMULATIONS 


115 


4.4 Growth of the Mean Flow 

Wygnanski, Champagne, & Marasli [58] (hereafter referred to as WCM) studied a 
set of small defect turbulent plane wakes created with a variety of wake generators. 
The generators were carefully designed such that the momentum thickness, 6$, was 
constant for all of the experiments. The Reynolds number for the various wakes they 
studied, based on freestream velocity and momentum thickness, ranged between 640 
and 3220 (which, assuming a Gaussian mean initial profile, corresponds to Reynolds 
numbers based on wake halfwidth and centerline velocity of between Re b = 589 and 
Re b = 2962, the high end of the computational range studied here.) They found that 
the far wake growth rates followed 



where x is the streamwise coordinate and x Q is a virtual origin. The coefficient A 0 
was found to vary depending on the particular wake generator used. The limiting 
values for A 0 that they observed were 0.270 for a solid strip set perpendicular to the 
flow direction, and 0.382 for a flat plate with a trailing edge flap which was oscillated 
at the frequency of the Karman mode. 

In order to compare with these experiments, straight dotted lines indicating the 
upper and lower bounds for the growth rates observed by WCM will be plotted 
in the following figures. These lines have been transformed into the appropriate 
computational variables using x — Xo = Uo{t — to)- Note that these lines have not 
been shifted to account for the virtual origins for the various flows. Only the slopes 
of the lines are significant. 

4.4.1 Effect of Disturbance Wavelength 

Figure 4.16 shows comparisons of the square of the wake halfwidths, b and 6, , versus 
non-dimensional time, t, for computations with a Reynolds number of Re b = 346 
and various combinations of initial disturbances. All the disturbance phases are zero. 
Up to a time of approximately t = 70, the growth rate of all the wakes shown is 
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(a) 



(b) 


Figure 4.16: Square of normalized wake halfwidth versus normalized time for Re = 
346 wakes with various combinations of initial disturbance wavelengths, o: 346(60)£— . 
O: 346(60)£° x . □: 346(60)°-. A: 346(60)** V : 346(60)**. (a) Direct halfwidthl. 
(b) Integral halfwidth b t . 
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dominated by the initial two-dimensional development of the Karman vortex street. 
The presence of disturbances other than the two-dimensional fundamental has little 

impact. 

Between a time of t = 70 and t = 200 the growth rate becomes highly dependent 
on the particular choice of initial disturbance. The most significant factor is the 
existence of a two-dimensional subharmonic disturbance. The wakes with the two- 
dimensional subharmonic (346(60fc 346(60)$*, and 346(60)$°) undergo a period 
of rapid spreading around a time of t = 100 which corresponds to the pairing of 
the spanwise rollers, while the wakes with only the two-dimensional fundamental 
disturbance grow more slowly. 

The late time growth rate, after a time of approximately t = 200 is determined by 
the wavelengths of the oblique disturbances. The wakes initiated with a subharmonic 
oblique disturbance, which develop strong three-dimensional motions, maintain a 
growth rate similar to the upper range observed by WCM. The wakes initiated with 
only a fundamental wavelength oblique disturbance, which are dominated by two- 
dimensional dynamics, have a late time growth rate well below the experimental!} 
observed level. The late time growth of these wakes is primarily due to the viscous 
diffusion of the coherent vortex street. 

4.4.2 Effect of Disturbance Phasing 

Figure 4.17 shows plots of the square of the wake halfwidths versus non-dimensional 
time for a set of three-dimensional simulation runs started from the 346(60)o r > * 1 dis- 
turbance condition (a two-dimensional fundamental and a three-dimensional oblique 
fundamental at 60°) overlayed on the results from the corresponding two-dimensional 
simulation (346(0)5”). It is readily apparent that the growth of these wakes is dom- 
inated by two-dimensional dynamics. This is not surprising given the fact that, as 
discussed previously, this choice of initial disturbances does not lead to development 
of strong three-dimensionality. As a result, the phase of the oblique disim nance has 
no significant effect on the growth of the mean flow. The late time growth rates for 
these wakes fall well below the levels observed by WCM. 
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(a) 



(b) 

Figure 4.17: Square of normalized wake halfwidth versus normalized time for 

346(60)oir X wakes (two-dimensional fundamental plus three-dimensional fundamen- 
tal at various phases). o: 346(60 ){£*. O: 346(60)^//. □: 346(60)^. : Cor- 

responding two-dimensional wake, 346(0)5". ( a ) Direct halfwidth b. (b) Integral 
halfwidth b,. 
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( a ) 



(b) 


Figure 4.18: Square of normalized wake halfwidth versus normalized time for 

346(60)oox XX wakes (two-dimensional fundamental and two-dimensional subharmonic 
plus three-dimensional fundamental at various phases), o: 346(60)gg*. 346(60) O o r ■ 

□ : 346(60)^. X . : Corresponding two-dimensional wake, 346(0 )5o*. (a) Direct 

halfwidth b. (b) Integral halfwidth b t . 
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(a) 



(b) 

Figure 4.19: Square of normalized wake halfwidth versus normalized time for 

346(60)o'^ x wakes (two-dimensional fundamental and two-dimensional subharmonic 
at ^ plus three-dimensional fundamental at various phases), o: 346(60)02*.. O: 

— jrj- — XT ^ 

346(60)01^. □: 346(60) o 2 i. x . : Corresponding two-dimensional wake, 346(0)5i x x . 

(a) Direct halfwidth b. (b) Integral halfwidth b t . 
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(a) 



(b) 

Figure 4.20: Square of normalized wake halfwidth versus normalized time for 

346(60)o ? i^ c wakes (two-dimensional fundamental and two-dimensional subharmonic 
at | plus three-dimensional fundamental at various phases), o: 346(60)o| T r - O: 

346(60)^. □: 346(60)J^. : Corresponding two-dimensional wake, 346(0)S|*. 

(a) Direct halfwidth b. (b) Integral halfwidth b,. 
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Figures 4.18, 4.19, and 4.20 show plots of the square of the wake halfwidths versus 
non-dimensional time for a set of three-dimensional simulation runs started from the 
346(60)o'^ disturbance condition (a two-dimensional fundamental, a two-dimension- 
al subharmonic, and a three-dimensional oblique fundamental at 60°). Here again 
the flow is dynamically two-dimensional and therefore the growth of the mean flow 
is very insensitive to the phasing of the three-dimensional fundamental disturbance. 
The phasing of the two-dimensional subharmonic determines the development of the 
mean flow almost entirely. The small deviations from the corresponding two-dimen- 
sional wakes at late times are due to the (relatively weak) three-dimensional structures 
which develop during the pairing of the spanwise rollers. However, these structures 
never dominate the flow dynamics. Again the late time growth of these flows is well 
below the range observed by WCM. 

Figure 4.21 shows plots of the square of the wake halfwidths versus non-dimension- 
al time for a set of three-dimensional simulation runs started from the 346(60)o^ r dis- 
turbance condition (a two-dimensional fundamental and a three-dimensional oblique 
subharmonic at 60°). In contrast to the wakes with oblique disturbances at the fun- 
damental wavelength, the late time growth of these wakes is substantially different 
from the corresponding two-dimensional case. The initial growth is still dominated 
by two-dimensional dynamics, but by a time of t = 125 three-dimensional processes 
clearly begin to take over, as evidenced by the substantial deviation from the two- 
dimensional growth curve. 

The effect of the phasing of the three-dimensional subharmonic on either halfwidth 
measure is minimal up to a time of t = 200. After a time of / = 200, however, the 
“direct” halfwidth measure shows what appears to be a strong phase dependence, the 
apparent dependence on phase is due to the fact that the mean flow profiles for these 
wakes are highly non-Gaussian, with several local maxima. The oblique disturbance 
phase has an impact on the details of the shape of the mean profile, as well as an 
effect on the maximum mean velocity. This has a large impact on the calculation 
of the “direct” halfwidth, b, even though the vorticity in the flows has spread over 
a similar extent. This is illustrated by the mean Ui velocity profiles shown in figure 
4.22. which correspond to the o and o lines in figure 4.21 at a time of approximately 
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(a) 



(b) 

Figure 4.21: Square of normalized wake halfwidth versus normalized time for 

346(60)oii I wakes (two-dimensional fundamental plus three-dimensional subharmonic 
at various phases), o: 346(60)o°J- C>: 346(60 ) O j T . □: 346(60) o J x . : Cor- 

responding two-dimensional wake, 346(0 )q**. (a) Direct halfwidth 6. (b) Integral 
halfwidth b,. 
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x2 


Figure 4.22: Mean velocity profiles for 346(60 )q xx and 346(60)o xx ;r wakes. : 

692(60)5 XX , t = 288.2. : 692(60)o rx r , t = 292.1. Arrows indicate locations for 

measurement of “direct” halfwidths. 

t = 300. The relative small differences in profile height and shape lead to a large 
difference in “direct” halfwidth b (see arrows). 

The actual effect of disturbance phase is relatively small as can be concluded 
from the plot of the integral halfwidth b t . What is most significant is that the late 
time growth rates for all of these wakes are similar and within the range observed by 
WCM. The presence of a three-dimensional subhaxmonic allows growth that is in line 
with natural wakes. This does not occur in the wakes that have a three-dimensional 
fundamental only. 

Figures 4.23, 4.24, and 4.25 show the square of halfwidths versus time for the 
346(60) o ^ / x wakes. Again, the oblique disturbance phasing has only a small impact 
on the growth of the mean flow. Though the magnitude of the late time width varies 
somewhat with the phasing of the oblique disturbances, the growth rate is not signifi- 
cantly impacted. Note also that so long as a three-dimensional subharmonic is present, 
the late time growth rate of the three-dimensional wakes varies little with the phase of 
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Figure 4.23: Square of normalized wake halfwidth versus normalized time for 

346(60)oox X wakes (two-dimensional fundamental and two-dimensional subharmonic 
plus three-dimensional subharmonic at various phases). o: 346(60 )g°x- 

346(60)oj r I - □: 346(60)£jL*. : Corresponding two-dimensional wake, 346(0)gj£. 

(a) Direct halfwidth b. (b) Integral halfwidth 6,. 
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(b) 


Figure 4.24: Square of normalized wake halfwidth versus normalized time for 

346(60)oi'] z wakes (two-dimensional fundamental and two-dimensional subharmonic 
at | plus three-dimensional subharmonic at various phases), o: 346(60 )oi x x - O: 

346(60)0^. □: 346(60) o | r . : Corresponding two-dimensional wake, 346(0)oi I r . 

(a) Direct halfwidth b. (b) Integral halfwidth b,. 
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(b) 

Figure 4.25: Square of normalized wake halfwidth versus normalized time for 

346(60)oi ? ^ x wakes (two-dimensional fundamental and two-dimensional subharmonic 
at - plus three-dimensional subharmonic at various phases), o: 346(60)oi r r . O: 

346(60)q|^. □: 346(60)qI x . : Corresponding two-dimensional wake, 346(0)g| x I . 

(a) Direct halfwidth b. (b) Integral halfwidth b t . 
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the two-dimensional subharmonic disturbance, or for that matter, with the presence 
or absence of the two-dimensional subharmonic (compare with the 346(60)o*x 1 wakes 
in figure 4.21*). The late time growth rates are again similar, though in some cases 
slightly above, the bounds observed by WCM. 

4.4.3 Effect of Reynolds Number 

Figure 4.26 show's a comparison of the square of the wake “direct” halfwidth versus 
time for a set of computations all initiated with the (?)(60 )q°* disturbance condition, 
but at different Reynolds numbers. Figure 4.27 are similar plots of the integral 
halfwidth. At early times (figure 4.26b and 4.27b), the low Reynolds number wake 
rolls up more slowly and peaks in width later than the other wakes. This is due to 
viscous effects. The 346(60)£°£ wake shows a similar, though much less pronounced 
effect. 

From the time the wake rolls up into the initial vortex street until a time of 
approximately t = 150, the wake growth is insensitive to Reynolds number. During 
this period the three-dimensional structures are growing in strength but have not yet 
become the dynamically dominant features in the wakes. After a time of t = 150, the 
growth begins to show some sensitivity to Reynolds number. The three-dimensional 
structures in the lowest Reynolds number wake, 119(60)^, are quickly diffused by 
viscosity and and therefore never become strong enough to overwhelm the two-dim- 
ensional dynamics of the early flow\ At higher Reynolds numbers, however, coherent 
three-dimensional structures dominate the late time growth. Thought there is some 
variation in the late time widths of the wakes as the Reynolds number is varied, 
the late time growth rate for all of the higher Reynolds number flows is similar. In 
fact the growth rates of the highest Reynolds number flows are somewhat lower than 
the growth rate for the 346(60)g°* wake. The reason for this is that the 346(60)g°* 
wake develops very strong three-dimensional structures that remain coherent for long 
periods. Those organized structures can efficiently work to spread the wake. At higher 
Reynolds numbers the coherent structures are disrupted, though not destroyed, by 
disorganized fine scale motions. The resulting reduced organization of the coherent 
structures slightly reduces the wake growth rate. 



CHAPTER 4. THREE-DIMENSIONAL SIMULATIONS 


129 



t 


(b) 


Figure 4.26: Square of normalized wake direct halfwidth versus normalized time for 
(?)(60)o°£ wakes (two-dimensional fundamental plus three-dimensional subharmonic 
with various Reynolds numbers). □: 119(60)o°*. o: 346(60)5°*. A: 692(60)5°*. V : 
1384(60)5°*. x: 2768(60)5°5- (b) is a magnification of the early time region of (a). 
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(a) 



(b) 


Figure 4.27: Square of normalized wake integral halfwidth versus normalized time for 
(?)( 60 )o°i wakes (two-dimensional fundamental plus three-dimensional subharmonic 
with various Reynolds numbers). □: 119(60)^. o: 346(60)^. A: 692(60)o°£- V : 
1384(60)ox?- x: 2768(60)5 °x- (h) is a magnification of the early time region of (a). 
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Figure 4.28: Comparison of halfwidth versus downstream distance for 119(60)5°* wake 
and experiments of Corke, Krull, & Ghassemi [14] (CKG). o: CKG experimental data. 

CKG curve fit. : 119(60)g£. 

4.5 Comparison with Experiment 

The 119(60)5°^ and 346(60)5°^ wakes in figures 4.26 and 4.27 has slight undulations 
in the late time growth curve visible in plots of both width measures. This starts 
at a time of approximately 80 with the peak which corresponds to the rollup and 
saturation of the Harman mode, and is followed by a period of neutral or negative 
growth, followed by growth again. This cycle of stronger growth followed by weaker 
growth continues as both wakes develop. 

Corke, Krull, &; Ghassemi [14] (hereafter referred to as CKG) studied this phe- 
nomenon in a Re = 119 spatially developing wake produced by a symmetric airfoil 
at zero angle of attack. They forced their wake with a two-dimensional disturbance 
at the fundamental frequency and a pair of oblique modes at the subharmonic in a 
manner similar to the computations presented here. They concluded that the fluc- 
tuations in the growth of the wake were due to a parametric resonance between the 
fundamental mode and the oblique subharmonic mode. The two-dimensional mode 
saturates and begins to feed energy into the three-dimensional oblique disturbance 
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through a secondary instability mechanism. The oblique mode in turn saturates, and 
the two-dimensional disturbance begins to grow again. 

Figure 4.28 shows a comparison between the direct wake width for the 119(60)5^ 
wake and data from figure 4 in CKG. The computed flow has been renormalized to 
match momentum thickness and shifted to put the virtual origin for the computed flow 
at x = 0. The coordinates were transformed using x/&ckg = Uo[t— to) /0 computation and 
( b/0)cKG = {b 1 6) computation- This rescaling does not account for differences in initial 
forcing magnitude between the experiment and the computation, or for the virtual 
origin of the experimental wake. Nevertheless, the match between the computed flow 
and the experimental data is quite reasonable. 


4.6 Selected Spectra, Mean Profiles, 
and Turbulence Statistics 

Figure 4.29 shows rescaled streamwise energy spectra (k x E(k x )) for the velocity 
fields in two of the computed wakes. The data has been plotted in this way in order 
to bring out regions of k x 5/ " 3 spectra which are expected in turbulent flows. In these 
coordinates, regions of k x 5 ^ 3 spectra will appear as horizontal lines. The spanwise 
and cross-stream (£(£ 2 ) and £(/c 3 )) spectra are very similar to the streamwise spectra 
and are not shown here. 

Figure 4.29a shows the rescaled streamwise energy spectrum for the moderate 
Reynolds number 346(60 )q°£ wake which appears in figure 4.13b. There is clearly 
no significant inertial (k x 7 , horizontal line in the plot) range in this flow. This is 
consistent with the fact that the flow is dominated by large coherent structures and 
has few small scales. This wake can not be considered turbulent. 

Figure 4.29b shows the streamwise energy spectrum for the highest Reynolds 
number wake computed, 2768(60)g°£. This is the wake which appears in figure 4.13b. 
This flow has a spectrum which is much more consistent with the expected spectrum 
of a turbulent flow. It shows a range in the spectrum that is close to the expected k x 5 ^ 3 
law (horizontal line). Since the range of ki over which the spectrum goes like k x 5 ^ 3 is 
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(a) 



(b) 


Figure 4.29: Rescaled streamwise energy spectra for 346(60)^ and 2768(60)g°* 

wakes. : (a) 346(60)g* t = 102.7. (b) 2768(60)g£, t = 102.4. 

Arrows indicate minimum grid wavenumber. 
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fairly short, the wake can not be considered fully turbulent. A proper classification 
for this wake is “transitional” . 

This “transitional” classification is borne out by the mean streamwise (ui) velocity 
profiles which appear in figure 4.30. (The symmetry of the initial conditions is such 
that the spanwise and cross-stream velocities have zero mean). These profiles corre- 
spond to the wakes which appear in figure 4.13. Only for the lowest Reynolds number 
case (119(60)5^ in figure 4.30a) where no significant three-dimensionality develops, 
is the mean profile Gaussian. For all of the higher Reynolds number wakes (figures 
4.30b-e) the coherent structures in the flow distort the shape of the mean profiles. 
Specifically, the “lumps” in the profiles that appear in the vicinity of x 2 = ±5.0 are a 
result of the strong streamwise structures that appear in figures 4.13b-e. The shape of 
the profile for the 2768(60)j££ wake (figure 4.30e), which is significantly non-Gaussian 
again indicates that that flow is not fully turbulent, but instead is still transitional. 

Figure 4.31 shows the second order velocity correlations, u'u'*. versus the cross- 
stream coordinate x 2 corresponding to the mean profiles in figure 4.30. Again, because 
of the symmetry of the initial conditions, u\u' 3 and u 2 u 3 are zero, so they are not 
shown. These profiles are typical of plane wakes, with the most intense unsteady 
motions occurring away from the wake centerline. As expected, the intensity of the 
unsteady motions increases with increasing Reynolds number. Once again the effects 
of the coherent structures are on the profiles in the vicinity of x 2 = ±5.0. This is 
particularly apparent in the u' 2 u' 2 profiles for the wakes in figures 4.31c-e. This implies 
that the coherent structures shown in figure 4.13 are still responsible for significant 
unsteadiness in the wakes. 

In a physical flow, the effects of coherent structures on long time averages are 
smeared out. This is due to the presence of low amplitude long wavelength motions 
that cause each set of coherent structures to be offset in space slightly from the other 
sets. Since the longest wavelength available in any of the simulations presented is 
not very large, the effective averaging domain is not large, and the effects of coherent 
structures on the averages are more pronounced than would be seen in a physical 
experiment, figures 4.32 and 4.33 are an attempt to get around this shortcoming. 
This has been done by taking mean data from several times near the time of interest 
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(approximately t i O.lf in this case) for each simulation. The separate mean data 
sets are normalized using the integral halfwidth at each time, 6 ,(t), and forming a 
combined mean. This effectively increases the sample size, and introduces the sort of 
jitter in the coherent motions that one would expect in a physical flow. 

The long time mean streamwise velocity profiles in figure 4.32 are considerably 
smoother than the single time mean profiles for the same datasets in figure 4.30. 
The profiles are significantly more Gaussian in appearance, though the effects of the 
coherent three-dimensional structures are still evident near the edges of the wakes. 

The same smoothing effect can be seen in the long time second order correlations 
in figure 4.33. The general shapes remain unchanged from the single time plots in 
figure 4.31, but the bimodal form of the ujuj, u' 2 u' 2 , and u 3 u 3 correlations are more 
clearly evident, particularly for the higher Reynolds number cases. 

Figure 4.34 shows single time mean streamwise velocity profiles for the same set of 
simulations shown in figure 4.30 but at a later time of approximately t = 200 for each 
wake. Again the single time means for all of the higher Reynolds number wakes sim- 
ulated are quite non-Gaussian due to the strong coherent streamwise structures that 
appear in the flow. Taking a long time average (figure 4.36) makes the mean velocity 
profiles appear more Gaussian, but the highest Reynolds number wakes continue to 
have distinctly non-Gaussian mean profiles. 

Figure 4.35 shows the single time second order velocity correlations for the same set 
of wake simulations, again at a time of approximately t = 200 for each wake. Figure 
4.37 shows the corresponding long time correlations. The most significant thing to 
note is that the Re = 119 wake (figure 4.35a) has significantly weaker unsteady 
motions than the wakes at Reynolds numbers. This is due to the fact that at that 
low Reynolds number, the wake never develops strong fine scale motions. 

The effect of fine scale motions is also apparent in the relative magnitudes of the 
second order velocity auto-correlations for the higher Reynolds number wakes. At 
the earlier time of t ~ 100 showm in figures 4.31 and 4.33 the it 2 u 2 correlation was 
approximately twice the magnitude of the ujtij and U 3 U 3 correlations. At t ~ 200, 
the onset of significant small scale motions have caused these three quantities to be 
both much larger in magnitude relative to the mean flow (note difference in scales in 
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figures 4.33 and 4.37) and much closer in magnitude to each other. This is due to 
the fact that at the earlier time, most of the velocity fluctuations are due to large 
scale coherent structures, which favor the u'^u'^ correlation. At the later time more of 
the velocity fluctuations are due to the more randomly distributed fine scales, which 
contribute equally to all three second order auto-correlations. 



Chapter 5 


Three-dimensional 
Topological Description 


In order to study the fine scale (high wavenumber), high gradient motions in the 
computed wakes, a topological classification method has been applied. The method 
is based on concisely summarizing the local flow structures in the space of the in- 
variants of the velocity gradient, strain rate, and rotation rate tensors. This allots 
information about the local flow geometry for every point in the incompressible flow 
to be presented in the form of two-dimensional joint probability density functions 
(pdfs) of various combinations of the invariants. These pdf plots facilitate the stud\ 
of global tends in the local structure of the fine scales in the velocity field. 

Only a small subset of the invariant space plots for the wakes that have been 
simulated will be presented here. See appendix E for a more complete set of plots. 


5.1 Topological Method 

An abbreviated description of the topological method used will be presented in this 
section. An extended description can be found in appendix D. 

The velocity gradient tensor 

Ajk = u h k (5-1) 
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can be split into a symmetric and an antisymmetric part 


Ajk — Sjk 4 " 14 j k 


(5.2) 


where the symmetric part 

Sjk = T 

is the rate-of-strain tensor, and the anti-symmetric part 


H j* = 7 }( u j,k ~ Ukj) 


(5.3) 


(5.4) 


is the rate-of-rotation tensor. The eigenvalues A of A 3 k satisfy the characteristic 
equation 

A 3 + PX 2 + QX + R = 0 (5.5) 

where P , Q, and R are referred to as the invariants of the tensor A : k and are given 
by 


P = -A k k = -Skk (5.6) 

<? = \{P 7 - An,A ki ) 

= - Sn.Su - W ik W tl ) ( 5 . 7 ) 

R = — detfA,*] 

= '-(-P^ + SPQ-AjkAuAtj) 

= l -(-P 3 + 3 PQ - S jk SkiS h - m) k W k ,S h ). (5.8) 

The values of the three invariants P, Q. and R completely determine the eigenval- 
ues. and therefore the eigenvectors, of the velocity gradient tensor A } k- Since the 
eigenvectors of Ajk determine the local flow kinematics, the local flow geometry is 
determined to within an arbitrary rotation by the location of the three invariants in 
( P.Q.R ) space. See Chong. Perry, and Cantwell [12] for a detailed discussion of the 
local flow geometries associated with the various regions of ( P,Q.R ) space. 
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Figure 5.1: Invariant Space for Incompressible Flows 


For incompressible flow 



P = 0 

(5.9) 

and the second and third i 

invariants reduce to 


Q 

= -\{S ]k S k} + W }k W kj ) 

(5.10) 

R 

= -^{SjkSkiSij + dWj k WkiSi 3 ). 

(5.11) 


Hence the local flow geometry is completely determined by the location of the second 
and third invariants of the velocity gradient tensor in ( Q,R ) space. 

The curve in {Q,R) space that separates characteristic equations with all real 
solutions (strain dominated local flow) from those with one real and two complex 
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solutions (rotation dominated local flow) is given by 

27i? 2 + 4g 3 = 0. (5.12) 

This curve, along with the Q axis, separates the ( Q,R ) space into four regions as 
shown in figure 5.1: 

• Above the separator and to the left of the Q axis, the local flow spirals in 
towards the local origin in a plane and then flow out along the third direction. 
This local flow geometry is referred to as a stable- vortex/stretching. 

• Above the separator and to the right of the Q axis, the local flow is toward the 
local origin along one axis and spirals out in a plane. This local flow geometry 
is referred to as an unstable- vortex/contracting. 

• Below the separator and to the left of the Q axis, the local flow approaches 
the origin along two axes and flows outward along the third. This local flow 
geometry is referred to as a stable-node/saddle/saddle. 

• Below the separator and to the right of the Q axis, the local flow approaches 
the origin along one axis and flows outward along the other two. This local flow 
geometry is referred to as an unstable-node/saddle/saddle. 

The flow geometry for a single data point can be determined simply by calculating 
the invariants and determining where they lie in ( Q , R) space. Global trends for large 
quantities of data can be examined by calculating the invariants for each point in the 
data set and constructing joint probability density functions (pdf's) in (Q.R) space 
for the entire data set. Moreover, motions with high gradients will tend to lie far 
from the origin in {Q,R) space, so the most intense motions will tend to separate 
themselves visually when viewed in the space of invariants. 

Other quantities related to the invariants of the velocity gradients tensor are also 
of interest. The second invariant. Q. may be written as the sum of quantities 

<? = -5(^4 + 

= Qs + Qu 


(5.13) 
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where 

Q, = (514) 

is the negative-definite second invariant of the rate-of-strain tensor and 

0 , = \w jk W, t ( 5 -' 5 ) 

is the positive-definite second invariant of the rate-of- rotation tensor. 

The second invariant of the rate-of-strain tensor is directly proportional to the 

mechanical dissipation rate 

9 = 2 vS ik S k j = -AuQ.. ( 5 - 16 ) 

Points with a large negative value for Q s have large dissipation. The second invariant 
of the rate-of-rotation tensor is equal to the enstrophy density 

MM = Q*. ( 5 - 17 ) 

Points with large positive values of Q w have high associated enstrophy density. Thus, 
the second invariant of the velocity-gradient tensor can be thought of as a measure of 
the relative importance of strain and rotation. Plots of — Q s vs. Q w reveal correlations 
between strain and rotation fields. Plots of Q s vs. R„ the third invariant of the rate- 
of-strain tensor 

R s = -\s jk S kl S h ( 5 - 18 ) 

reveal trends in the type of rate-of-strain field associated with high dissipation re- 
gions. Since the rate-of-strain tensor is symmetric it only has real eigenvalues, hence 
its second and third invariant will always fall underneath the separator, and the pos- 
sible rate-of-strain topologies will be limited to either stable-node/saddle/saddle or 
unst able-node/saddle/saddle. 
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Finally, the vortex stretching rate, cr, may be expressed in terms of the invariants 
of the velocity-gradient and rate-of-strain tensors 

o = W }k W kl S h = R s - R. (5.19) 

Plots of the enstrophy density versus stretching. Q w vs. a can be examined to reveal 
global trends. 

5.2 Effect of Initial Conditions 

Figure 5.2 shows plots of the joint pdf of Q and R for a set of Ret, = 346 wakes with 
a variety of combinations of initial disturbance wavelengths. The first feature that is 
apparent is that, as one would expect from the discussion of the wakes in physical 
space presented earlier, the wake initiated with only a two-dimensional fundamental 
and a three-dimensional fundamental (346(60)°^ in figure 5.2a) has fine scale motions 
which have much lower gradients (are much weaker) than the other wakes (recall that 
high gradient motions tend to appear far from the origin in invariant space). Further, 
though the wake which has an additional two-dimensional subharmonic disturbance 
(346(60 )°q^ in figure 5.2e) has gradients that are on par with the wakes that have 
a three-dimensional subharmonic disturbance (figures 5.2b-d), those high gradients 
account for a much smaller percentage of the flow as evidenced by the relative small 
area inside the second contour level. Note that the contour levels shown are logarith- 
mic. with each contour level being ten times the value of the previous level, thus large 
differences in the area inside of a given contour level equate to very large differences 
in the relative volume of the fluid in physical space represented. 

The second feature that is readily apparent is that the overall shape of the joint pdf 
of Q and R is the same for all of the wakes, with high gradients motions favoring the 
upper left (stable- vortex/stretching) and lower right (unstable- node/saddle/saddle) 
regions of (Q.R) space. This same “skewed teardrop” shape has been observed by 
Soria ef al. [42] and Chen et al. [11] in studies of the gradients in compressible and 
incompressible temporally evolving mixing layers, Blackburn, Mansour. k Cantwell 





CHAPTER 5. THREE-DIMENSIONAL TOPOLOGICAL DESCRIPTION 161 



00 

o 


I- 


b0 

05 

C 

1— 1 



"oj 

II 

> 


ft 

**o 

£ 


0) 

ss 

u 

a 

© o 

H O 

ft 

o 


o 

U ( 


P 

o 



oc 

CO 

"6 

„ , 

_ . 

o 

ft 




*C 




tn 

.2 

00 

*c 

It 

(t 

II 

> 

■w 

i— ! 




i 

H H 
o o 

CO 

H O 

o 

o 


;c 

ft 


£ 





CO 



CO 

c 

II 

' — ' 

-c 


QJ 

QC 

QS 



o 

c 

C^J 

t-HH 


GJ 

II 





> 


o> 

H H 
O H 
HO 

t*— 


0 

o 

U-* 


-c 


Cu 

-o 


CO 

c 




c 

JD 





CHAPTER 5. THREE-DIMENSIONAL TOPOLOGICAL DESCRIPTION 162 


[3] in a computation of a turbulent channel flow, and Sondergaard et a 1. [41] in 

a survey of computations of compressible and incompressible shear flows including 
mixing layers, wakes, and isotropic and sheared homogeneous turbulence. 

It is apparent from the wide variety of flow fields which exhibit this feature that 
it is a characteristic of solutions to the Navier-Stokes equations (the pdf of Q and 
R generated from random gradient fields do not exhibit this feature). Cantwell [6] 
performed an analysis of the evolution of the velocity gradient tensor in incompressible 
flows and found that, given certain assumptions, the second and third invariants were 
restricted to an attractor in invariant space with a shape similar to the one which has 
been observed. 

Figure 5.3 shows plots of the joint probability density functions of the second 
and third invariants of the rate-of-strain tensor ( Q s and R s ) for the same wakes 
shown in figure 5.2. Once again, the much weaker gradients in the two-dimensional 
fundamental / three-dimensional oblique fundamental case (346(60 )q^ in figure 5.3a) 
are apparent. And again, the 346(60)°5J wake (figure 5.3e) shows gradients on par 
with the simulations which have a three-dimensional subharmonic, but with a much 
smaller volume of high gradient motions occurring. 

All of the pdf s have the same general shape, with highly dissipating motions 
(large negative Q s ) showing a strong preference for the unstable-node/saddle/saddle 
type topology. This strong preference continues even for more moderately dissipating 
regions as demonstrated by figure 5.4, which is a magnified version of figure 5.3b. 
These moderately dissipating motions account for the majority of the total dissipation 
in the flow. Again, these same trends were observed by Chen et ah, Soria et al.. 
Blackburn et al.. and Sondergaard et al. in a wide range of other flows. 

Figure 5.5 shows plots of the joint pdf of -Q s (dissipation) and Q w (enstrophy 
density) for the same set of simulations shown in figures 5.2 and 5.3. Once again the 
much lower gradients for the three-dimensional fundamental cases (figures 5.5a and 
5.5e) is apparent. These wakes also have a pdf shape which which is different than 
those for the wakes with a three-dimensional subharmonic (figures 5.5b-d). The “L" 
shaped pdf in figures 5.5a and 5.5e, which tends to hug the axes, is common in non- 
turbulent flows and flows in which the turbulent motions have not had time to fully 





CHAPTER 5. THREE-DIMENSIONAL TOPOLOGICAL DESCRIPTION 165 



Figure 5.3: Contour plots of joint pdf of Q, vs. R s for Ret, = 346 wakes with various initial disturbance wavelengths, 
(a) 346(60)g”, t = 196.9. (b) 346(60)52*. t = 204.8. (c) 346(60)55}*. t = 181.1. (d) 346(60)^°. t = 197.0. 
(e) 340(60)55”. t = 202.8. 
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0.0 R s 



Figure 5.4: Enlarged plot of pdf of Q s vs. R s for 346(60 )q°*. t = 204.8. Corresponds 
to figure 5.3b. 

develop. The highest mechanical dissipation rates in such flows tend to occur in the 
high strain-rate regions between large organized vortical structures, thus the shape of 
the pdf. The half ellipsoid shape seen in figures 5.5b-d is much more typical of st rongly 
three-dimensional flows, w'here regions of high strain-rate are closely intermingled with 
regions of high vorticity. 

The “tongues" of highly rotational points that appear in the flows with a three- 
dimensional subharmonic, particularly in figures 5.5b and 5.5c, correspond to the 
centers of the strong streamwise structures seen in physical space (since those struc- 
tures contain the highest enstrophy density motions in the wake). It is clear from 
figure 5.5 that while the most dissipating motions (largest |(?,|) tend to be associated 
with the regions of highest enstrophy density (scattered points near —Q s = Q w line, 
most obvious in figures 5.5a.b and e), the bulk of the high enstrophy density points in 
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-0.025 0.0 0.025 

Figure 5.6: Contour plots of joint pdf of Q w vs. R s — R for 346(60)qoo • t = 197.0. 

the wake have only moderate dissipation rates. This leads to the conclusion that the 
bulk of the most mechanically dissipating motions in the wake occur in the vicinity 
of the highest enstrophy density motions, but are in fact associated with more mod- 
erate enstrophy density levels. This is not entirely unexpected, since many models of 
the structure of turbulent motions produce peak dissipation rates separated from the 
peak enstrophy density location (c.f. Burger's vortex). 

Figure 5.6 shows the joint pdf of Q u , (enstrophy density) and R s — R (vortex 
stretching rate) for the 346(60 )qqq wake. The points with the highest enstrophy 
density occur in regions of the flow that have the strongest vortex stretching. This 
comes as no surprise since the effect of stretching is to amplify the vorticity aligned 
with the strain field. This result is typical of all the three-dimensional flows which 
have been studied. 
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Figure* 5.7: Sop caption page 172. 
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dimensional subharmonic with various Reynolds numbers), (a) Rc h = 346, t = 204.8. (b) Re b - 692, t = 201.9. 
(r) Re b = 1384, t = 199.3. (d) Re b = 2768, t = 194.6. 
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5.3 Effect of Reynolds Number 

Figure 5.7 illustrates the effect of Reynolds number on the invariant pdf s. Shown is 
the joint pdf of Q and R for (?)(60)q°x wakes all at approximately the same develop- 
mental time. The effects of increasing the Reynolds number is to greatly increase both 
the magnitude of the highest gradients and the fraction of the flow which contains 
high gradient motions. 

Changing the Reynolds number has minimal effect on the overall shape of the 
pdf’s (given that the Reynolds number is high enough to allow for strong three-dim- 
ensional motions). The only effect of increasing the Reynolds number is to make the 
highest gradient regions appear a bit more scattered (compare figure 5.7d and figure 
5.2b on page 159). 

5.4 Time Evolution in Invariant Space 

Figure 5.8 illustrates the time evolution of the wake flows in invariant space. Shown 
is the joint pdf of Q s and R s for the 346(60)££ wake. The initial shape of the pdf 
is determined by the initial conditions of the flow (figure 5.8a), but it rapidly relaxes 
to the general shape which has been observed in the other flows, and which will be 
maintained throughout the remainder of the evolution (figure 5.8b). 

The strong gradient motions, which are associated with strong three-dimension- 
alitv. become much more significant (note the much larger area encompassed b\ the 
second and third contour levels in figure 5.8c as compared to figure 5.8b). The 
magnitude of the highest gradient motions then begin drop as the wake approaches a 
self similar regime where all the gradients should decay like At the time shown 

these wakes have not yet entered the self similar decay regime. 

This developmental cycle is common to all of the invariants for all of the three- 
dimensional wakes simulated, and is independent of Reynolds number. 
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ntermediate strain-rate. 
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5.5 Rate- of- Strain Distribution and 

Vorticity- Strain Alignment 

Figures 5.9 and 5.10 show plots of the probability density functions of the normal- 
ized principal rates-of-strain, a, 0, and 7 , for the 346(60)^ and 1384(60)5°* wakes 
respectively. The strain rates are sorted such that they are in descending order 

o > 0 > 7 (5.20) 

and the rate-of-strain in each principal strain direction has been normalized by the 
magnitude of the intermediate rate-of-strain, \0\. Thus the normalized intermediate 
rate-of-strain. 0 , can take on values of positive or negative one only 

0 = ±1.0. (5.21 ) 

This is the reason for the delta function distribution for 0. 

There is a direct correspondence between the sign of 0 and the local rate-of-strain 
topology. If 3 — +1.0. then the local rate-of-strain topology is of the type unstable- 
node/saddle/saddle. If 0 = —1.0, then the local rate-of-strain topology is of the type 
stable-node/saddle/saddle. 

If the pdfs are formed using all of the grid points in the computational domain 
(including those in the freestream), they appear as in figures 5 . 9 a and 5 . 10 a. Approx- 
imately two-thirds of the points have a rate-of-strain field consisting of two positive 
and one negative rate-of-strain (0 = 1.0). The other two rates-of-strain have broad 
distributions with peaks around a ss 2.0 and 7 ss — 3 . 0 . 

If the pdf's are conditioned so as to include only those points with high mechanical 
dissipation, which is proportional to — Q s , then they appear as in figures 5 . 9 b and 
5.10b. The high dissipation points have rate-of-strain fields which almost exclusively 
consist of two positive and one negative rate-of-strain (0 = 1 . 0 , unstable-node/saddle/ 
saddle rate-of-strain topology as was apparent from figure 5 . 3 b). The distributions of 
q and 7 still have rather broad distributions, but with more distinct peaks at a % 1.5 
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and 7 ~ -2.5 This result is insensitive to Reynolds number, and holds for all of the 
three-dimensional wakes computed, even those with relative weak three-dimensional 
motions (e.g. 346(60)o*J ). 

Ashurst et a L [2] studied a direct numerical simulation of isotropic turbulence 
and found that the strain-rates for the most dissipating motions were in the ratio of 
a : 3 :') = 3:1: - 4 , which they speculated might be a universal ratio. The studies of 
Soria et a 1. [42] and Chen et al. [11] and Sondergaard et a I [41] found that the ratio 
of principal rates-of-strain depended on the specific flow examined. The strain-rate 
ratio of a : l3 : 7 1.5 : 1 : - 2.5 found here for the three-dimensional incompressible 

plane wake adds weight to the conclusion that the ratio is indeed flow-dependent 

Figures 5.11 and 5.12 show plots of the probability density function of cosine of 
the angle between the vorticity vector and each of the three principal strain-rate 
directions for the 346(60)^ and 1384(60)^ wakes. The pdf of the cosine is plotted 
because in those coordinates randomly distributed three-dimensional vorticity vectors 
will result in uniform pdf's. 

There is a clear tendency for the vorticity vector to be aligned (cosine equals 
one) with the intermediate rate-of-strain (/?) direction and nonaligned (cosine equals 
zero) with the most compressive rate-of-strain ( 7 ) direction in both wakes (figures 
5 . 11 a and 5 . 12 a). When only the points with the highest dissipation rates, and hence 
the most intense local rate-of-strain fields, are included in the pdf, this tendency is 
strongly enhanced (figures 5 . 11 b and 5.12b). The vorticity aligns almost exclusively 
with the intermediate rate-of-strain direction and is nearly always approximately 
perpendicular to the most compressive rate-of-strain direction. 

Again, similar results were observed by Ashurst et a 1. [2] for isotropic turbulence. 
They also found that the vorticity tended to align with the intermediate rate-of-strain 
direction. The studies of Soria et al. [42] and Chen et al. [11], Sondergaard et al. 
[41], Blackburn et al. [3] and Tsinober et al. [47] found the same tendency for a 
wide variety of other flows. The evidence that this is a universal characteristic of 
turbulence is becoming very convincing. 
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Conclusions 


6.1 Numerical Method 

The numerical method developed to perform the simulations used for this study has 
proven to be both effective and efficient. For planar flows with compact vorticity 
fields, matching of the velocities at the edge of a finite sized computational domain 
to irrotational flow solutions is a workable alternative to nonlinear mapping of the 
infinite direction to a finite domain. 

The velocity matching method provides uniform resolution of the vortical region of 
the flow without wasting large fractions of the available grid points on the freestream 
or over-resolution of the centerplane of the flow. Use of a Fourier method on the 
uniform grid also allows for changes in resolution to be accomplished by simply trun- 
cating the transformed data or padding it with zeros at runtime. No special routine 
or extra computational time is needed to rebuild the grid or dataset after a resolution 
change. 

The addition of a growing grid in the cross-stream direction allows the computa- 
tional domain to adapt to the changing size of the vorticity field as the simulation 
evolves. This permits the flow to remain resolved even over long simulation times 
with minimal need for intervention on the part of the user. Changes in grid size of a 
factor of two or more while the simulated wake grew in extent by factors of five were 
routine when running the simulations for this study. 
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6.2 The Incompressible Plane Wake 

6.2.1 Physical Development 

The results presented here indicate that oblique disturbances at the subharmonic 
wavelength are very important to the development of strong three-dimensionality 
in the temporally evolving incompressible plane wake. Simulations started with a 
two-dimensional fundamental disturbance and a three-dimensional disturbance at the 
subharmonic wavelength developed strong streamwise structures stretching between 
corrugated spanwise rollers on the same side of the wake. These structures create 
very high rates-of-strain which lead to the appearance of fine scale motions and rapid 
growth rates in the far wake. Simulations started with three dimensional disturbances 
at only the fundamental wavelength remain almost two-dimensional, with relatively 
weak rate-of-strain fields, no significant small scale motions, and sluggish late time 
growth patterns which nearly match the corresponding purely two-dimensional wakes. 

Varying the Reynolds number affects the intensity and scale of the structures 
in the flow. This effect is especially strong in wakes which develop strong coherent 
three-dimensional structures (wakes initiated with an oblique subharmonic). At high 
Reynolds number these wakes develop very intense fine scale motions even at moderate 
Reynolds numbers. The most intense fine scales tend to appear in the vicinity of the 
coherent three-dimensional structures. Wakes which do not develop strong three-dim- 
ensional coherent structures (wakes initiated with an oblique fundamental) develop 
small scales only at the very highest Reynolds numbers. The Reynolds number of 
the wake has minimal impact on the growth rate of the far wake. W 7 akes with very 
low Reynolds number do not develop strong three-dimensional structures, and grow 
somewhat more slowly at late times. Wakes with moderate to high Reynolds numbers 
grow at similar rates given similar initial conditions. 

The addition of longer wavelength disturbances allows for scale changes in the flow 
and can temporarily augment the growth rate of the wake. This has implications for 
the interpretation of temporal simulations, where the flow is restricted to a maximum 
wavelength corresponding to the longest wavelength in the initial disturbance field. 
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At late times in the simulation the flow will not follow a growth path which can be 
related to that of a spatially developing wake. 

This result also has implications for experiments involving spatially developing 
wakes. Any experiment is by necessity limited by the dimensions of the facility 
in which it is performed. The sensitivity of the wake to subharmonic disturbances 
which has been observed in these computations suggests that the rate of growth of 
experimentally studied wakes may be sensitive to the low frequency spectral content 
of naturally occurring disturbances in experimental facilities. 

6-2-2 Topological Development 

Topological analysis of the fine scale, high gradient motions in the incompressible 
wakes revealed that the wakes with a three-dimensional subharmonic have a greater 
quantity of more intense high gradient motions as compared to wakes with a three- 
dimensional fundamental only. Changes in the phase of the initial disturbances were 
found to have minimal impact on the overall distribution of the high gradient motions. 
Increasing the wake Reynolds number increased the intensity of the gradients while 
preserving the shape of the pdf’s in invariant space. 

All of the three-dimensional wakes simulated had joint probability density func- 
tions which were similar to those observed in other three-dimensional flows: 

• Joint pdf's of the second and third invariants of the velocity gradient tensor have 
a characteristic "skewed teardrop” shape, with high gradient motions tending 
to be of stable- vortex /stretching and unstable-node/saddle/saddle topological 
types. 


• Joint pdf’s of the second and third invariants of the rate-of-strain tensor indicate 
that the most dissipative motions are associated exclusively with an unstable- 
node/saddle/saddle type strain topology. More moderately dissipating motions, 
which account for the majority of the dissipation in the flow, are also very 
strongly associated with an unstable-node/saddle/saddle type strain topology. 
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• Joint pdfs of the enstrophy density and vortex stretching indicate that the 
highly rotational motions occur in regions where the vortex stretching is posi- 
tive, even at late times when the intensity of all gradients are decreasing in the 
wake. 

Examination of the rate-of-strain distributions and the vorticity-strain alignment 
indicates that highly dissipating motions have a strong tendency to have two positive 
and one negative principal rates-of-strain (which corresponds to an unstable-node/ 
saddle/saddle type rate-of-strain topology) in the ratio of 1.5 : 1 : -2.5. These 
motions also tend to have the vorticity vector aligned with the intermediate (positi\e) 
strain rate direction and nearly perpendicular to the most compressive (negative) 
strain rate direction. These results appear to be insensitive to both Reynolds number 
and initial condition. 


6.3 Future Work 

As a direct extension of this study, the following work is recommended: 

• The effect of the angle of the oblique disturbance needs to be examined, at least 
to the extent of verifying the broad range of amplified secondary instabilities 
predicted bv the work of Flemming [16]. 

• With the recent availability of larger and faster parallel machines, the existing 
code should be used to simulate wakes at higher Reynolds numbers and for 
longer times. This is necessary to verify that the profound difference between 
wakes with short wavelength oblique disturbances and wakes with long wave- 
length oblique disturbances continue to exist at very high Reynolds numbers 
and at late times. 

• To adequately simulate late time behavior, the wake simulations need to be 
run with the addition of longer wavelength disturbances to allow for continued 
unconfined growth. This should also allow for much realistic calculation of late 
time mean profiles and turbulence statistics. 
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In addition, the code which was developed for this study can readily be used to 
study a wide variety of planar free shear flows, including mixing layers (as demon- 
strated in section 2.8.1), skewed mixing layers, and momentumless wakes. Applica- 
tion of the code to these flows requires only that the proper initial conditions for the 
vorticity field and free stream flow velocities be defined. 



Appendix A 

Classical Similarity Theory 


This appendix will present a brief review of the arguments and conclusions of sim- 
ilarity theory as applied to free shear flows in general and the incompressible plane 
wake in particular. 


A.l Preliminaries 

A. 1.1 Energy Transport Equation 

The Navier-Stokes equations for incompressible flow may be written as 


Ujj — 0 


Uj,t + (ujUl + ~Sji - =0 

where S 3 k is the rate-of-strain tensor defined by 


Sjk = 9 + Ufc o) • 


(A.l) 
(A. 2) 


(A. 3) 
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Multiplying equation A. 2 by u k gives 


u k 


U h t + ^ Sji — 


= {ujUk) t + ^UjUkUi -f ^ ( u k 6ji + Uj6 k i) — 2u[u k Sji + 

- 2 Zs jk + 2v (Sj, (S kl + Wki) + S kl (S jt + Wji)) = 0 
P 

where W jk is the rate-of-rotation tensor defined by 


(A.4) 


Wjk — 0 { u j,k u k,j , 


(A. 5) 


Taking the trace of equation A.4 and noting 5„ = 0 for incompressible flow gives 
the transport equation for the kinetic energy E = u k u k /2 


(A. 6) 


+ £14. + - u k - 2 ^Uj 5 j> + $ = 0 

P J.k 


where <J> is the dissipation of (total) kinetic energy 


$> = 2vS ik S ki > 0. 


(A. 7) 


The flow variables may be split into mean and fluctuating parts 

u j = Uj+u' j p = p + p (A. 8) 

where the overbar signifies an appropriate average of the variable (temporal or en- 
semble for spatially developing flows, spatial or ensemble for temporally developing 
flows) and the primed quantities are deviations from that mean. 

Using this decomposition, and applying the averaging procedure to the Navier- 
Stokes equations A. 2 yields the Reynolds equations 


(A. 9) 
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iij't + (ujUk + -Sjk — 2 vSjk 4 — 0 (A. 10) 

V P P ) m 

with a mean kinetic energy defined by 

E = (ujUj + u jU} ) /2 

= £ + ? (A-ll) 

where E is the energy of the mean flow and q 2 is the turbulent kinetic energy. 

A transport equation for the energy of the mean flow may be derived from equation 
A. 10 by multiplying by u : and rearranging 

E t + ^IJkE + — n* — 2uHjSjk + u j u 'j u 'k^j A -vSjkSkj ~ u j u k ^ A ~ ^ (A. 12) 

Similarly, a transport equation for the total mean energy can be derived from 
equation A. 6. 

+ <? 2 ,, + ( Eu k + q 2 Uk + ~Tk - 2vu ] S lh + ZZju'u'*. 4- q 2 u' k + ^u' k 

- 2»C^Tk) h + 2vSjkSk] + 2«/S£s£ = 0. (A.13) 

Subtracting out the equation for the mean flow' energy, equation A. 12, from the 
above equation yields an equation for the transport of mean turbulent kinetic energy 

+ feu* + 7< + ~ p < - 2u^S])j = 11 “ ^ (A ‘ 14) 

where 

n = -u' 3 u' k s jk (A. 15) 

is known as the turbulence production, and 

v = 2v'5£5i~ 


(A. 16) 
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is the dissipation of turbulent kinetic energy. The production, 11 serves to transfer 
kinetic energy from the mean flow to the unsteady flow, while the dissipation serves 
to convert turbulent kinetic energy into heat. 

It has been observed in a number of flows with simple geometries that at suffi- 
ciently high Reynolds number the production scales locally with the dissipation 

If ^ ip. (A. 17) 


or 

-^ k S ]k ~2vS^ (A. 18) 

A complete theoretical understanding of this relationship between production and 
dissipation is still lacking. The evidence for equation A. 18 is empirical and confined 
to a limited range of flows. However this is used as the basis for arguments presented 
in the following sections. 

A. 2 Scales of Motion 

Equation A. 18 allows for an estimate of the relative size of the large and small scales 
in the flow. Assume that the velocity fluctuations and the mean velocity both scale 
with the same reference velocity 


u'j ~ Uj ~ U 0 (A. 19) 

and that the largest scales in the flow scale with some reference length 6. This allows 
an estimate of the scale of the left hand side of equation A. 18 

n = ~^7 k s ]k ~ L -j- (A. 20) 

The scaling of the right hand side is more problematic, since it involves the mean 
of products of the gradient of the fluctuations. Clearly the appropriate velocity scale 
to use is l o since u ' ^ Uo* but the appropriate length scale to use is unclear. 
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There are two turbulence length scales, or microscales, that are commonly used to 
complete the above scaling: The Taylor microscale, and the Kolmogorov microscale. 

A. 2.1 The Taylor Microscale 

The Taylor microscale, A, is defined such that the proportionality given in equation 
A. 18 is satisfied 

(a. 2i) 

6 A 2 

giving 

This microscale represents an upper bound on on the range of scales that contribute 
to significant turbulent kinetic energy dissipation. 

In this study, two version of the Taylor microscale were calculated for use in 
estimating the resolution of the simulations that were performed. The first version is 
one that appears in many experimental studies where obtaining the three-dimensional 
velocity and velocity gradient field is very difficult. It relies on the measurement of 
only one velocity (usually u^) 

(A. 23) 

The second version used here includes the velocities and gradients in all three 
directions 

which is derived from the definition which appears above for the general case. 

The Taylor microscale also arises in correlation functions which appear in the 
theory of isotropic turbulence. For a complete description see the discussions in 
Hinze [18]. 
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A. 2. 2 The Kolmogorov Microscale 

The second scale is a result of introducing both a new turbulence length scale and a 
new turbulence velocity scale so as satisfy the proportionality in equation A. 18 


t / 0 3 v 2 

^ 1 / 


(A. 25) 


To close the definition, the Reynolds number of the resulting scale is chosen to be one 



v 


(A. 26) 


These lead to the relations 

(A. 27) 

and 

(A. 28) 

The Kolmogorov microscale represents a lower bound on the on the range of scales 
that contribute to significant turbulent kinetic energy dissipation. For both the Taylor 
and Kolmogorov microscales it should be kept in mind that they are only estimates 
of the actual turbulent length scales. 




A. 3 Evolution of the Spatially Evolving 
Incompressible Plane Wake 

A. 3.1 Momentum Balance 

Figure A.l shows schematically the wake behind a symmetric body with a drag force 
per unit span of D . 

For a sufficiently large control volume (the dotted box in figure A.l), only the 
downstream edge of the control volume and the drag on the body contribute to 
the momentum balance in the control volume. The upstream boundary and top 
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3*2 = H / 2 



x 2 = -H/2 - 

Figure A.l: Schematic of spatially evolving wake. Control volume for momentum 
balance. 

and bottom of the control volume are at free stream conditions, and therefore ha\c 
negligible contribution to the momentum. Thus, for sufficiently large H, an integral 
momentum balance over the control volume becomes 


D 

P 


W/2 r 


/ 

-Hj 2 


Ul(3 , 2)(Wl(J2 


Hoc) T 


'p(x 2 ) - Po 


Til (2*2) ~ T lloc 

P 


dx 2 


(A. 29) 


where r tJ is the stress tensor. 
In the far wake. 


^CC 

Combining this with the similarity assumptions 


Uqc ~ H 1 (3*2 ) _ ^ f £2 \ 

U 0 ~ U \6 ) 


(A. 30) 


(A. 31) 


and 


( P(£ 2 ) 


Til (X 2 ) - Tile 




P 


P 


(A. 32) 
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x 2 — H / 2 



Figure A. 2: Schematic of spatially evolving wake. Control volume for energy balance. 


where S is a measure of the mean flow width and U 0 = u* — u<£ is a reference velocity, 
allows the momentum balance to be written as 


D 

/™X 


n I to n J 2d 

/ " ( t ) ( t ) + u ° 6 J f 


-H/26 


<T 


(A. 33) 


In the far wake, the second integral is much smaller than the first., and the momentum 
balance can be approximated as 



(A. 34) 


l nder the assumption of self-similarity (the profile is a function of x 2 /S only), the 
integral is a constant, in this case of order one. Thus, in the far w r ake 



(A. 35) 
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A. 3. 2 Mean Kinetic Energy 

The evolution of the mean kinetic energy can be derived from the integral of equation 
A. 12 over the infinitesimal control volume shown in figure A. 2 

J [£, ( ] dV + J [ (u k E + ^u k - 2 i/UjSjk + + 2 j vSjkSkj]dV 

= -Jpy k S jk \dV. (A.36) 

V 

Assuming the flow is stationary, the first integral is zero (constant total energy in 
the control volume). The second integral may be converted into a contour integral, 
giving 

J [tJ fc £ + -u k - 2v\ijS jk + Ujit'jU'f. + 2vS jk S k jjdVt k 
n P 

= - J [ujt i' k Sj*]dV. (A-37) 

v 

where fl is the contour surface and Q k is a surface normal. 

Let AF = I'oc ~ «7 and Ap = p ■<* - P ■ Evaluating over the volume in A. 2 and 
assuming the mean flow is in the Xj direction yields 



= - j\u',u' t S, k \dV. 

V 


(A. 38) 



APPENDIX A. CLASSICAL SIMILARITY THEORY 


196 


Introducing the similarity assumptions 


A U . (x 2 

U 0 ~ U \ 6 


(A. 39) 


A p _ 3A U 
P ~ 2 






= h (*2 


(A. 40) 


( A .4 1 ) 


(A. 42) 


Po , 3 r 


■ 


/ /(?)-(? 


-H/26 


H fit n/zc 

+ m, / / KfMf) (A - 43) 


-Hf 25 


The first term on the left is zero since, by equation A. 35, f/o<5 is constant. The third 
term dies off much more rapidly in the far wake than either of the remaining terms, 
so it may be dropped. This leaves the approximate relation 


n ( zc n / zo 

«-(«*)., / / K t)" 


(A. 44; 


-H/26 


The integrals are constants of order one, hence 


Coo (CoO, ~ 'o- 


(A. 45) 
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Combining this with equation A. 35 gives 




u'J 2 UW 


and 



the asymptotic growth laws for the spatially evolving incompressible plane 


A. 4 Evolution of the Temporally Evolving 
Incompressible Plane Wake 

A. 4.1 Momentum Balance 

For the temporally evolving plane wake the total momentum deficit per 
area M (momentum per unit span per unit length) is a constant. Thus 


P 


H/ 2 

= I (u<x - U\{x 2 ))dx 2 . 
-HI 2 


Applying the similarity assumption from equation A. 31 yields 

H/ 26 

= U 0 

P 


H ( 26 

M .. f .(X2\,(X 2 \ 


-H/26 


Again, the integral is a constant of order one, thus 


M 

U 0 6 ~ . 

P 


(A. 46) 


(A. 47) 
wake. 


unit plan 


(A. 48) 


(A. 49) 


(A. 50) 
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A. 4. 2 Mean Kinetic Energy 

The development of the results for the mean kinetic energy are identical to those for 
the spatially evolving wake except that the time derivative of the mean energy is not 
dropped. The equivalent to equation A. 38 is 



= -/[«>! S )t pV. (A. 51) 

V' 

Applying the similarity assumptions from equations A. 39 through A. 42 gives 



-HI 21 -H/ 26 


The first term is zero since H is a constant. The same holds for the second term 
since, by equation A. 50, L oS is a constant. All of the terms which are differentiated 
with respect to X\ are zero since in a temporal wake there is no spatial variation of 



APPENDIX A. CLASSICAL SIMILARITY THEORY 


199 


the mean. This leaves the third term on the left and the term on the right 

H/26 H/26 

^ir I * 2 (f M?) = / KfMf)' 

-H/26 -H/26 


Again, the integrals are constants of order one, hence 

(USD, 


2 —«• 


Combining this with the result from equation A. 50 yields 


and 


s ~ j 

p i 


Uo ~ 

\4 p) 


(A. 53) 


(A. 54) 


(A. 55) 


(A. 56) 


the asymptotic growth laws for the temporally evolving incompressible plane wake. 





Appendix B 

Linear Stability Theory 


In the interest of completeness, this appendix will present a brief overview of linear 
stability theory. A brief outline of the methodology used in the present study to 
generate the disturbance eigenfunctions will also be presented. 


B.l Mean Flow 

The mean base wake flow for all of the present simulation runs was the parallel 
Gaussian profile (properly non-dimensionalized) 

€/, = 1 - (B.l) 

ii 2 = u 3 = 0 (B.2) 

p = const (B.3) 

where Au<£o is the dimensionless centerline velocity deficit, and C) is a scaling factor. 
The scaling factor C] was chosen to be 0.69315 and Au<j*o was chosen to be 0.692 

which gives an initial wake halfwidth 6 0 of 1.0, and an initial Reynolds number based 

on halfwidth of 0.692///. These values were used in the previous experiments of Sato 
Kuriki [39]. Corke, Krull. Ghassemi [14], and in the computations of Chen. 
Cantwell. Mansour [10]. 
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B.2 Linearized Disturbance Equations 

Starting with the full incompressible, uniform density Navier-Stokes equations 


Co 

II 

O 

(B.4) 

p.j 

Uj't + u k u Jt k 4- — = vuj'kk 

(B.5) 

VM — puk,iui,k 

(B.6) 

a total flow consisting of the mean base flow defined in equations 

B.l through B.3 

plus a small perturbation 

v.j — U\8\j -|- u'j 

(B.7) 

p = p + p' 

(B.8) 

is substituted into equations B.4 through B.6. 

= 0 

(B.9) 

P j “h p f ' 

+ Uj ft + + U k ){u\,kf>lj + U jk ) + — 

+ «U) ( B1 °) 

P,kk + P[ k k = + U' ]k )(ui t j6ik + u' k j) 

(B.ll) 

The mean flow satisfies the Navier-Stokes equations, hence the 

mean flow terms 

may be subtracted out from equations B.9 through B.ll leaving the nonlinear distur- 

bance equations 

<,■ = 0 

(B.12) 

p f 

v! ] t + ! + u' 2 Ui t 2 ^\] + u 'k u '],k + = vu 'j ,kk 

(B.13) 

P,kk = Pftu 1,2 u 2,1 "1" U k,l U l,k) 

(B.14) 
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The perturbations are assumed to be small and products of perturbation terms are 


dropped, giving the linearized disturbance equations 

ujj = 0 (B.15) 

P f 

u ',,t + j 4- u' k ui <k 6ij + — = vu ' kk (B.16) 

J J p 

p[kk = 2/»Si,2«2,i (B.17) 

The disturbance flow is assumed to have the form of a traveling wave 

u] = u J (x 2 )e , (° I1+/Jl3 - c<) (B.18) 

p' = p(x 2 )e ,(oil+/3l3 ~ cf) (B.19) 


where q and 0 are real wavenumbers which determine the wavelength and wave angle, 
and c = c# + ic t is a complex wavespeed. 

Substituting into the linearized disturbance equations 


laiti + U2,2 + i0u$ = 0 (B.20) 

Z Q 

i(aui - c)ui + ii 1 , 2 u 2 = p + ^(uj, 22 - (q 2 + 0 2 )ui) (B.21 ) 

P 

Z ( Q U i - c)u 2 = — ~P,2 + ^(^ 2,22 ~ (« 2 + Z? 2 )^) (B.22) 

P 

z(qU] - c)u 3 = p + z / (u 3,22 - (o 2 + ^ 2 )u 3 ) (B.23) 

p 

P.22 - (a 2 + P 7 )p = 2zopux, 2 u 2 (B.24) 


Equations B.20 through B.23 may be combined to eliminate the pressure. After 
some manipulations, equations B.20 through B.23 can be reduced to the equations 

v ( D 2 — (a 2 + /? 2 )) — i(ahi — c ) (.D 2 — (a 2 + /3 2 )) + ta(D 2 u\) u 2 = 0 (B.25) 
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and 

v D 2 — (a 2 + /? 2 )) — — c)] (au 3 — f)u\) = — /?(Z)t2i)u 2 . (B.26) 

Here the differentiation operator D = (),2 has been introduced for the sake of clarity. 

Any physically meaningful disturbance must decay to zero as |x 2 | — ► oo, thus the 
appropriate boundary conditions for these equations are 

[u 2 , £>ii 2 ] 0 as |x 2 | — ► oo (B.27) 

and 

oit 3 - 0ui — ♦ 0 as |x 2 | — ► oo (B.28) 

Equation B.25 is the well known Orr-Sommerfeld equation. Equation B.26 is 
known as the Squire equation. The form of equations B.25 and B.26 indicate that 
there are are tw r o classes of solutions to the set of linear equations B.20 through 
B.23 ([8]). The first class is comprised of solutions to equations B.25 and B.27 with 
equations B.26 and B.28 used only to solve for tii and u 3 once tt 2 is known. The 
second class is comprised of solutions with u 2 = 0 which satisfy equations B.26 and 
B.28. The solutions we will seek here are in the first class. 


B.3 Solution of the Linear Equations 

The technique used to solve equation B.25 subject to B.27 was the spectral method 
developed by Spalart et a i. A brief overview of the procedure will be given here. 
The reader is directed to Spalart, Moser, and Rogers [44] for a complete detailed 
description. 

Spalart et a J. begin by defining a vorticity perturbation component perpendicular 


to the wave vector (a, (3) 



, — + 0^3 

UJ — — 

(B.29) 


y/a 2 + P 2 

where 

u>! = Du3 — 

(B.30) 
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and 

u> 3 = iau 2 — Du\ (B31) 

are the perturbation vorticities. 

Applying continuity (equation B.20) allows u? 1 to be written in terms of u 2 only 
i\ja 2 + fPd' 1 = \D 2 — (a 2 + /3 2 )] (B.32) 

Substituting equation B.32 into equation B.25 gives 

\v [D 2 - (q 2 + 8 2 )) - i(qui - c)] \Ja 2 + = —a(D 2 ui)ti 2 (B.33) 

with boundary conditions 

u? 1 — ► 0 as \x 2 \ — ► oo. (B.34) 

Equations B.33 and B.34, along with B.32, are solved using a standard Galerkin 
method. The known and unknown functions, /(x 2 ), are expanded in terms of a set of 
orthogonal basis functions R 3 derived from the (1,1) Jacobi polynomials on a mapped 


coordinate r/ 


u x {rj) = ajRjirj) 

(B.35) 

u 2 {rj) = 6jv_iF/v_i(77) + VrV(7/) + bjRj(ri) 

(B.36) 

Uliv) = CjRjil) 

(B.37) 

Rj(v) = (1 - ri 2 )P?' 1 ] (v) 

(B.38) 

i = unh (!) 

(B.39) 

where x” is a scaling factor for the mapping and are exponential 

determined from 

“extra” functions 

[. D 2 - (a 2 + /J 2 )] r, = R, 

(B.40) 


which come about from the inversion of the Poisson equation B.32. 
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Derivatives of the basis functions Rj(ti) can be expressed in terms of the basis 
functions themselves using the recursion relation 

DR > - + + 

= L, k R k (B.41) 


Substituting equations B.35, B.36, and B.41 into equation B.32 allows the unknown 
coefficients a : to be written explicitly in terms of of the unknown coefficients bj 


vv + W{x° 2 y 

Mjkbk 


+ (ej — (a:2) 2 (o 2 + /? 2 )) b : 

+ /j+2&j+2, 

l dj-ibj-2 + ( x 2) 2 ^2’ 


j = 1, ...N-2; 
j = N-l, N; 

(B.42) 


where d,, e ; , and fj are known rational functions of the expansion index j. 

Substituting the above into equation B.33 and applying a scalar Galerkin method 
with test function Rk results in an N x TV matrix eigenvalue problem for the coefficients 
Oj of the form 

EjkCik = ca 3 (B.43) 

where E } h is a complicated matrix (which will not be presented explicitly here) formed 
from Ljk (the derivative recursion matrix), M } k (which relates the unknown coeffi- 
cients a } to of the unknown coefficients bj), the known expansion for u 1? and the 
known parameters a and (3. 

This eigenvalue problem is solved using a standard numerical package which re- 
turns the N complex eigenvalues, c, and their corresponding eigenvectors. The most 
unstable eigenfunction for the given input parameters (tij, a, and /?) is the one corre- 
sponding to the eigenvalue with the largest complex component C{. That eigenfunction 
is normalized and used to form the needed disturbance function. 



Appendix C 

Aliasing and Alias Control 


This appendix presents an overview of numerical aliasing and the techniques used in 
this study to control the associated errors. 


C.l Discrete Fourier Transforms and Aliasing 

The one-dimensional discrete Fourier transform (DFT) of a series a 3 of length N is 
given by 



(C.l) 

A T -1 

Oj = £ a n e 2 -^. 

n=0 

(C.2) 


with inverse transform 


Consider a series a } which represents a complex sinusoid with integer wavenumber 
k and constant amplitude C 

a : = Ce 2 *'». (C.3) 

The DFT of this series is 


A’-l 


a n = £ V = C6((k - n) mod N) 

A j=o 


(C.4) 
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where 6 is the discrete delta function 


6U) = 



if j =0 

otherwise. 


(C.5) 


Hence a wave with wavenumber k in physical space is transformed into Fourier space 
with a wavenumber ( k — n) mod N . For example, if k = TV L, the series is trans- 
formed as if it were a series with a wavenumber of n = L = k — N. This effect can 
been seen directly in physical space since for any integer M 


a, = Ce 2r '» = Ce *”> 2Li± F Un . 


(C.6) 


Thus data separated by M periods TV in physical space are the same. Wavenumbers 
k + MN are said to be “aliased” to wavenumber k. 

Aliasing errors occur in practice when operations on data with wavenumber span 
N increases the wavenumber span to greater than TV. Such is the case with the 
non-linear terms in the Navier-Stokes equations, which are bi-linear products of the 
dependent variables. The remainder of this discussion will be limited to such products. 

To examine aliasing in bi-linear terms, consider the one dimensional product of 
two series of length N 


N—1 N - 1 

V' i e 2 *'^ 


Cj = ajb } = a n t 2 *'K Y. b ” 


n—0 

N- 1 


m=0 


Y a n b m e 27r ' . 


n=0 m= 0 


(C.7) 


This product has an unaliased wavenumber span of 2TV — 1, nearly twice the span of 
a j and b r The length TV DFT of c 3 is 


c k = Y H a n b m 6{(n + m - k) mod TV). (C.8) 

n=0 m— 0 

Hence modes with (n + m) > N are aliased to modes with k = (n + m) — N. 
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The standard numerical DFT algorithms in use today have a span N which is 
even (typically a power of 2 with optionally one or more factors of 3) and span a 
wavenumber space of (1 — y) < & < y* Products of such series span a wavenum- 
ber range of (2 — N) < k < N. Figure C.l illustrates aliasing with such a modal 
arrangement. 

Multiplication 



Figure C.l: Example of aliasing 


In order to properly handle products of length N series, length N alias free prod- 
ucts must be formed. The alias error from such products can be eliminated in one of 
two ways, by truncation and by phase shifting. These methods are outlined below. 


C.2 Dealiasing using Truncation 

Since the product of two series each with a wavenumber span of N results in a series 
with a wavenumber span of (2N — 1 ) a simple method of obtaining a alias free product 
is to use transforms of length 2 N. The length N transforms a } and bj are padded 
with zeros to form length 2N transforms and transformed to physical space to form 
length 2N series. The series are used to form the needed products and transformed 
back to wave space with a length 2 N transform. An alias free length N transform 
can then be extracted by discarding modes beyond the desired span in wavespace. 

For the case of a wavenumber space (1 — y) < k < y identical results can be 
obtained using transforms of length rather than 27V , with resulting savings of time 
and storage space. This can been understood by considering that the worst cases for 
aliasing are when two modes of wavenumber y are combined, giving a product which 
is aliased to (A' — ~ N ) = — y and w T hen two modes of wavenumber (1 — y) are 
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combined, giving a product which is aliased to (|JV — N — 2) = (y — 2), both of 
which are outside of the desired alias free wavenumber space of(l — y)<A:<y. 

In practice, the maximum transform length N is determined by the limitations 
imposed by the available computer hardware, so the alias free wavenumber space 
carried is in practice 1 — y < & < y . Since two thirds of the available transform length 
N can be made alias free using truncation, this method of dealiasing is commonly 
referred to as the 2/3 rule. 


C.3 Dealiasing using Phase Shifts 

Equation C.8 can be rewritten in the form 


c : = X X a n i m . (C.9) 

n+m=j n+m=j±N 

> ' > v ' 

alias free alias error 

If the transforms are performed on a shifted mesh (which manifests in wavespace as a 
multiplication of each Fourier mode aj by a phase factor e' jA ) and then shifted back 
to the original mesh the results are 


* — IjA 

C J — e 


* inAt _tmA . \ ' * inAI J 

/ ^ a n t b m e + 2_j a n e 

\n+m~j n+m=j±A r 

= X a n b m + e ±lNA X 

n+m=j n + mrrj’i A r 


(C.10) 


alias free 


alias error 


The unaliased part of the product is unaffected by the phase shift while the alias 
error part is multiplied by a phase factor e ±tNA . This can be used to exactly eliminate 
the aliasing error by evaluating c 2 on two meshes shifted by half a cell width from 
each other. Then 

e ±tNA 2 _ e ±«N(A,±£) _ _ e ±»NA, (C.l 1 ) 


and the alias free product is simply the average of the two evaluations of c r 
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If exact dealiasing is not essential, as is the case for most time accurate simulations, 
multiple evaluations at each time substep can be avoided by evaluating the bi-linear 
products on a shifted grid at at each time substep. By choosing a random phase shift 
e ±iNA ever y other time substep, and using the half-cell offset shift e ±t/v ( A± n ) f G r 
the subsequent time substep aliasing error is eliminated to the same order as the time 
advance algorithm. This random shift method is used in the present code to control 
one-dimensional aliasing errors. 


C.4 Multidimensional Dealiasing 


Multi-dimensional Fourier transforms are obtained by applying separate one dimen- 
sional transforms in each direction 


*mnq. 


MNQ 


M - 1 N - 1 0-1 

LEE 

j =0 k = 0 1=0 




- 2 tti ( d 


■S+ ■) 


(C.12) 


The transform in each direction is independent of the other directions, hence aliasing 
in each direction may be treated independently. 

For the three-dimensional DFT, each pair of modes in a bi-linear product can 
combine to form one of four types of terms: an unaliased term, a term aliased in one 
direction, a term aliased in two of the three directions, or a term aliased in all three 
directions. Hence 


c £={j,U} “ 


+ 


+ 


^ " &v t i={m ,n,g} ^n= {r,a,t } 

rTi+flr: 

k 

s V y 

alias free 

) * dm “I" } ^ A } ^ hyftbfi 

k+{±M,0,0) *+{0,±N,0} *+{0,0,±Q} 


aliased in one direction only 

^ drn^n H" ^ ^ ^ni^n ^ ^ d m^n 

rfi -f + = 

fc+{±M,±N,0} *±{±M t 0,±Q} k±{0,±N.±Q} 

" 

aliased in two of three directions 
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+ 


rft+* = 

k±{±M,±N,±Q } 


aliased in all three directions 


(C.13) 


Dealiasing in multiple dimensions can be accomplished in the same ways as in 
one dimension, with truncation, shifting, or a combination of the two. Using only 
truncation is inconvenient since it reduces the number of useful modes to less than one 
third of the available modes. Dealiasing by pure phase shifting is also inconvenient 
since it would require two evaluations for each direction dealiased, for a total of six 
evaluations per time substep. 

The present code uses a combination of the two methods. Following Patterson 
and Orszag [31], modes with 



(C.14) 


are truncated. This eliminates the two and three dimensional aliasing as per the 
2/3 rule, leaving only the one dimensional aliasing. Truncating in this way increased 
the useful modes to approximately half of the total modes. Rogallo [34] advocates 
the somewhat less severe rule of truncating only those modes which are aliased in 
more than one direction (as opposed to the ellipsoid in C.14 which also truncates 
some modes which are aliased in only one direction) however the increased number 
of useful nodes (to approximately two thirds of the total number of nodes) is not 
significant when the resulting shape of the useful wavespace is considered. 

The remaining one dimensional aliasing errors can be handled using phase shifting 
as described above for the one dimensional case. The transform is evaluated twice on 
grids shifted by half a cell width in each direction, and the alias free product is the 
average of the two evaluations. 

Since exact dealiasing is not required for the present code, random phase shifting 
in each direction at alternate time substeps as outlined above for the one dimensional 
case is used to cancel out the aliasing error to the order of the time advance algorithm. 
The remaining aliasing error appears as a small amplitude random forcing function 
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at low wavenumbers [34]. So long as the computed flow remains well resolved (hence 
the energy in the highest wavenumber motions is several orders of magnitude below 
the most energetic wavenumbers) any residual aliasing error will be at least several 
orders of magnitude below the energy of the alias free solution. 



Appendix D 

Topological Classification 


Topological methods are useful in the description of vector fields and are coming in- 
creasingly into use as a means to study the large data sets produced by numerical 
simulations. Chong, Perry, and Cantwell [12] have carried out a general classification 
of the various types of linearized three-dimensional flows which can occur in com- 
pressible and incompressible flow'. This classification method was used by Cantwell, 
Chen, and Lewis [7] and Chen, Cantwell and Mansour [9] to analyze the the topology 
of flow structures in experimental measurements of a pulsed low-speed diffusion flame 
and direct numerical simulations of a compressible plane wake. Chen et. al. [11] and 
Soria et. al. [42] used this method to study the small scale motions in numerical 
simulations of a variety of compressible and incompressible flows including w'akes, 
mixing layers, isotropic turbulence, and homogeneous shear flows. 

The method is based on concisely summarizing local flow structures in the space 
of invariants of the velocity gradient tensor. In these studies, the velocity gradient 
tensor, A tJ = u tJ is calculated at each point in the flow, and the invariants of the 
velocity gradient tensor, as well as the invariants of the rate-of-strain and rate-of- 
rotation tensors are calculated. Plots of the joint probability density functions (pdf’s) 
of the invariants for the entire flow reveal global trends in the geometry of the velocity 
field which would be difficult if not impossible to discover using other techniques. 
They also allow' the study of how structures in invariant space (which correspond to 
specific local flow geometries) correspond to structures observed in physical space. 
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There are several good reasons for studying the velocity gradient tensor as op- 
posed to the raw velocity field. Primary among them is that results are coordinate 
independent (invariant under any affine transformation) and independent of the frame 
of reference of the observer. Moreover, in the case of incompressible flow where the 
first invariant is zero, the three dimensional physical field, which can be infinite in 
extent, can be represented in a finite region of the two dimensional space of the second 
and third invariants. Finally, large scale motions are associated with relatively low 
gradients while small scale motions are associated with high gradients. Thus different 
length scales in the physical flow naturally tend to be sorted into different regions in 
invariant space, and thus may be examined separately. 

While the present study focuses on the velocity field, it should be noted that this 
method may also be applied to any smooth vector field of interest. These can include 
the vorticity field, the scalar gradient field, or the pressure gradient field. 

D.l Local Flow Trajectories 

The instantaneous trajectory of any fluid particle in a flow field is determined by the 
solution of the convection equation 


Xj ,t = u j (D.l) 

where x 2 is the location of the given fluid particle and u 2 is the flow velocity evaluated 
at the particle position. 

For fluid in the neighborhood of some reference fluid particle at location x c } moving 
with the local flow velocity u c ., equation D.l may be used to obtain an equation for 
the local relative flow. 

Xj.t ~ x ),t = u j ~ u ] 

= {Xk - X c k )u 2ik + (x* - x*)(x, - X, C )Uj,fci + . . . 

~ U Jt k(x k - X k ). 


(D.2) 
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Hence, in a frame of reference moving with some particle in the flow, x'- = Xj — x c j, the 
trajectories of fluid particles near the reference particle with respect to the reference 
particle are determined by the solutions to the linear equation 

X 'j,t = U hk X 'k- ( D ' 3 ) 

The flow in this local frame of reference is completely determined by the eigen- 
vectors of the velocity gradient tensor 

Aj k = Uj.k (D.4) 


the symmetric part of which 

Sjk = + Uk, 3 ) (D.5) 

is the rate-of-strain tensor, and the anti-symmetric part of which 

Wjk = 2 ~ u kj) ( D - 6 ) 

is the rate-of-rotation tensor. If Ajk has only real eigenvalues, the local flow is strain 
dominated and the local flow consists of fluid moving inward to or outward from the 
origin along distinct axes.. If A } k has a pair of complex eigenvalues, the local flow is 
rotation dominated and the local flow consists of fluid spiraling around the origin in 
one plane and flowing inward or outward along the third direction. 


D.2 Eigenvalues 

The eigenvalues, A, and complex eigenvectors, e : , of Ajk satisfy the eigenvalue equa- 
tion 


\Ajk A<5 jfcjej — 0 


(D.7) 
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where the eigenvalues A are solutions to the characteristic equation 


det[A ifc - A S jk ] = 0. (D.8) 

Equation D.8 expands to the third order algebraic equation 

A 3 + /W 2 + QA + R = 0 (D.9) 

where P, Q , and R are referred to as the invariants of the tensor Aj k and are given 
by 


P = ~A kk = -S kk (D.10) 

Q = - A ]k A k] ) 

= - S, k S ki - W jk W ki ) (D. 11) 

R = — det[i4jjt] 

= } -(-P 3 + 3PQ-A jk A kl A h ) 

= l -(-P 3 + 3 PQ - SjkSuSu - 3 WjkWuSij). (D.12) 

The set of solutions for A from equation D.9 fall into one of three categories: All 
three A are real and distinct; All three A are real and at least two of them equal; Or 
one A is real and the other two are complex conjugates. The category which a given 
set of A falls into is completely determined by the location of the three invariants in 
( P , Q , R ) space. 

The surface in (P, Q , R) space that separates characteristic equations with all real 
solutions from those with one real and two complex solutions (and hence the surface 
on which the characteristic equation has three real solution with at least two equal) 
is given by 

27 R 2 + (4 P 3 - 18 PQ)R + (4 Q 3 - P 2 Q 2 ) = 0 (D.13) 
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A detailed discussion of the properties of this surface and a guide to the solutions 
for the resulting sets of (and hence the local flow geometry) that can occur in the 
various domains in ( P,Q,R ) space is given in Chong, Perry, and Cantwell [12]. 


D.3 Incompressible Flows 


For incompressible flows, 

P = Ajj = 0 

(D.14) 

and equations D.ll and D.12 simplify to 


Q = + W, t W k ,) 

(D.15) 

R = -LSjtSuSij + SW^WuStj). 

o 

(D.16) 


The surface which divides real from imaginary solutions described by equation D.13 
simplifies to the curve 

27i? 2 + 4Q 3 = 0. (D.17) 

Hence, for incompressible flows, the local flow geometry is completely determined by 
the location of the Q and R invariants in (£?,/?) space. 

Figure D.l shows all of the possible local flow geometries for incompressible flow. 
Below the curve given by equation D.17 all three eigenvalues of A } k are real, and the 
flow is dominated by strain type motions. To the left of the R = 0 axis, the local flow 
has fluid moving inwards toward the origin along two of the principal directions and 
outward from it along the third principal direction (stable node/saddle/saddle). To 
the right of the R = 0 axis, the local flow has fluid moving outward from reference 
point along two of the principal directions and inward towards the reference point 
along the third direction (unstable node/saddle/saddle). 

Above the curve given by equation D.17 two of the eigenvalues of A tJ are complex, 
and the local flow is dominated by rotational motions. To the left of the R = 0 axis, 
the local flow has fluid spiraling inwards toward the reference point in a plane and 
moving outward from it in the third direction (stable focus/stretching). To the right 
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of the R = 0 axis, the local flow has fluid moving toward reference point in one 
direction and spiraling outward in a plane (unstable focus/contracting). 

D.4 Joint Pdf’s of Invariants 

This relatively simple mapping of the location of the Q and R invariants in ( Q , R) 
space to the local flow geometry allows for a very concise summary of local geometry 
of a large numbers of points in a flow field. Q and R are calculated for each point of 
interest (often the whole flow) and the resulting data is presented as a joint probability 
density function (pdf) in the invariant plane. This allows the local flow geometries to 
be studied in a global framework which allows overall trends to be easily distinguished. 
This technique is especially useful for studying the smallest scales in the flow, since 
small scales will have high gradients and their invariants will tend to lie far from the 
origin in ( Q. R ) space where any trends can easily be distinguished. 

D.5 Enstrophy Density, Dissipation, 
and Vortex Stretching 

Other quantities related to the Q and R invariants are also of interest. In particular, 
Q and R maybe be split into contributions from the rate-of-strain and rate-of-rotation 
tensors. These contributions can be directly related to physical quantities in the flow. 

The second invariant of the velocity gradient tensor can be rewritten as the dif- 
ference of two pos-i-tive-def-inite quantities 

Q = -\(S»Sj t - 

= Qs + Qw (D.18) 

where Q s is the second invariant of the rate-of-strain tensor Sjk and Q w is the second 
invariant of the rate-of-rotation tensor W 3 k . This decomposition allows the the study 
of the relative importance of strain and rotation (enstrophy density) in the local flow' 
geometry. If Q is large and positive, then the local enstrophy density is large and 
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dominates the strain ( Q w ^ Qs)* If Q is large and negative, then the local strain 
is large and dominates the enstrophy density ( Q s ^ C?™). Plots of joint pdf of Q s 
versus Q w will reveal correlations between strain and rotation in the flow. 

The third invariant, R can be similarly split 

r = A(s jk s kl s h + w lk w kl s h ) 

= R, - a (D.19) 

where R, is the third invariant of the rate-of-strain tensor Sjk and a = WjkWkiSij 
represents stretching of vorticity. Plots of the joint pdf of u = R s — R and the other 
invariants will reveal correlations between vortex stretching and other quantities. Of 
particular interest is the correlation between vortex stretching and the enstrophy 
density (2 Q w ). 

The invariants of the rate-of-strain tensor Sjk are of additional interest. For in- 
compressible flow, the rate mechanical dissipation of energy due to viscosity is related 
to the second invariant of Sjk by 


9 = 2 vSjkSkj = -4 vQ 3 . (D.20) 

Hence strongly dissipating regions in a flow have large negative values of Q,. Since 
Sjk is symmetric, all its eigenvalues must be real, hence its invariants lie under the 
real-imaginary dividing curve given by equation D.17. By studying plots of Q s versus 
R s correlations between the mechanical dissipation and the local flow geometry can 
be examined. 
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(a) 346(60)g£, t = 196.9. (b) 346(60)g£. t = 204.8. (c) 346(60)5°*. t = 181.1. (d) 346(60)5°°. t = 197.0. 
(e) 346(60)^5*. < = 202.8. 
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Figure E.2: Contour plots of joint pdf of Q, vs. R , for Reb = 346 wakes with various initial disturbance wavelengths, 
(a) 346(60)°” t = 196.9. (b) 346(60)g£. t = 204.8. (c) 346(60)5°*. t = 181.1. (d) 346(60)5™ t = 197.0. 
(e) 346(60)^5*. t = 202.8. 
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wavelengths, (a) 346(60)°**, t = 196.9. (b) 346(60)^. t = 204.8. (c) 346(60)5°*. t = 181.1. (d) 346(60); 
/ = 197.0. (e) 346(60)°^. t = 202.8. 
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Figure E.4: See caption page 233. 
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Figure E.4: See caption page 233. 
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Figure E.4: Contour plots of joint pdf of Q w vs. R, - R for Re h = 346 wakes with various initial disturbance 
wavelengths, (a) 346(60)°^, t = 196.9. (b) 346(60)g£. t = 204.8. (c) 346(60)^. t = 181.1. (d) 346(60)5™ 
t = 197.0. (e) 346(60)°^. t = 202.8. 
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Figure E.8: See caption page 241. 
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Figure E.8: Contour plots of joint pdf of Q w vs. R, - R for <?)(60)S£ wakes (two-dimensional fundamental plus 
three-dimensional subharmonic with various Reynolds numbers), (a) Re b = 346, t = 204.8. (b) Re b = 692, 

t = 201.9. (c) Re b = 1384, t = 199.3. (d) Re b = 2768, t = 194.6. 
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Figure E.9: Time evolution of contour plots of joint pdf of Q vs. R for 346(60)5°* wake (two-dimensional funda 
mental plus three-dimensional subharmonic), (a) t = 22.8. (b) t = 52.8. (c) t = 102.7. (d) t = 204.8. 
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Figure E.12: See caption page 249. 
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ure E.12: Time evolution of contour plots of joint pdf of Q w vs. R, - R for 346(60)^ wake (two-dimensional 
damental plus three-dimensional subharmonic), (a) t = 22.8. (b) t = 52.8. (c) t = 102.7. (d) t — 204.8. 
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Figure E.15: Time evolution of contour plots of joint pdf of — Q a vs. Q w for 1384(60)5°* wake (two-dimensional 
fundamental plus three-dimensional subharmonic), (a) t = 28.4. (b) t, = 51.1. (c) t = 101.5. (d) t = 199.3. 
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Figure E.16: See caption page 257. 
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gure E.16: Time evolution of contour plots of joint pdf of Q w vs. R a - R for 1384(60)^ wake (two-dimensional 
ndamenta.1 plus three-dimensional subharmonic), (a) t = 28.4. (b) t = 51.1. (c) t = 101.5. (d) t = 199.3. 
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Figure F.2: See caption page 263. 
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Figure F.2: Mean velocity profiles and second order turbulence profiles for (?)(60)5°* wakes (two-dimensional 
fundamental plus three-dimensional subharmonic with various Reynolds numbers), (a) Re b = 346, t = 204.8. 
(b) Ile b = 692, t = 201.9. (c) Re b = 1384, t = 199.3. (d) Re b = 2768, t = 194.6. 
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Figure F.3: Mean velocity profiles and second order turbulence profiles for 346(60)5°* wake (two-dimensional 
fundamental plus three-dimensional subharmonic). (a) t — 22.8. (b) t = 52.8. (c) t = 102.7. (d) t = 204.8. 
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Figure F.4: See caption page 267. 
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Figure F.4: Mean velocity profiles and second order turbulence profiles for 1384(60)j£* wake (two-dimensional 
fundamental plus three-dimensional subharmonic), (a) t = 28.4. (b) t = 51.1. (c) t = 101.5. (d) t — 199.3. 







Appendix G 

Wakes with an 
Oblique Fundamental 


As discussed in the introduction, Williamson [52] and Williamson & Prasad [53, 54, 55] 
suggested an alternate to the mechanism examined in the bulk of this study for the 
development of strong three-dimensional motions in the plane wake. Their experi- 
ments indicated that strong, highly oblique coherent structures could be produced by 
an interaction between fundamental wavelength disturbances shed at a small oblique 
angle from the wake generator (as is common for the wakes behind bluff bodies such 
as cylinders) and long wavelength two-dimensional waves w'hich grow due to the hy- 
drodynamic instability of the far wake. 

In order to verify that this proposed mechanism is indeed a legitimate path for 
the development of three-dimensional structures, a small set of simulations were run 
to examine the behavior of such wakes. The simulations w-ere initiated with the 
same initial mean profile as used in the earlier simulations (see section 2.7 on page 
36). Fundamental disturbances at angles of between 5° and 30° with respect to 
the spanwise direction (0° being the non-oblique two-dimensional fundamental used 
in the earlier simulations) and two-dimensional (0°) subharmonic disturbances were 
superimposed on the mean. In each case, the width of the computational domain w'as 
chosen such that the fundamental disturbance was periodic in both the streamw’ise 
and spanwise directions. Note that because the fundamental disturbance is oblique. 
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there is no relevant measure of phase between the fundamental and the subharmonic. 
The phase, as defined for the two-dimensional wakes in chapter 3, effectively varies 
from 0 to 27 r along the span of the simulations. 

Figures G.l, G.2, and G.3 show oblique views of an iso-enstrophy density contour 
for wakes initiated with the oblique fundamental at 5°, 15 , and 30 respectively. The 
Reynolds number of these wakes, Rcb — 180, was chosen to approximately match the 
experiments of Williamson et a 1. The times (all approximately t ~ 100) match the 
times for the wakes in figure 4.13 in section 4.3.2 on page 111. The viewing angle is 
also the same as for the earlier simulations. 

The interaction of the oblique fundamental and the two-dimensional subharmonic 
causes the rollers of the (oblique) Karman vortex street which initially forms from 
the oblique fundamental disturbance to wrap around one another. This reorients 
some of the spanwise vorticity from the Karman vortex street into the streamwise 
direction. At low angles this wrapping is a result of the essentially two-dimensional 
interaction of the fundamental and the subharmonic. Each spanwise location sees a 
different phase between the fundamental and the subharmonic, thus each spanwise 
section pairs in a slightly different way as per the discussion in section 3.3.2. This is 
illustrated by the cuts through the vorticity field of the 5° wake which are shown in 
figure G.4. Eventually the rollers become sufficiently entwined for the three-dimen- 
sionality of the flow to become significant. At higher angles, the same basic process 
occurs, but the three-dimensionality of the flow begins to effect the dynamics at an 
earlier time. 

As Williamson et a 1. suggest, this mechanism is a likely candidate for the cell 
pattern seen in the far of the plane wake by Cimbala, Nagib, & Roshko [13]. Figure 
G.5 is a top view of the 30° oblique wake in figure G.3 (the free stream flow is from 
left to right). The cell pattern here roughly matches the pattern seen in figure 19 of 
Cimbala et a 1. 

The higher the angle of the oblique fundamental, the more quickly the rollers of 
the Karman vortex street become distorted, and the more quickly streamwise vorticity 
is produced. Only at the highest angle shown, 30 , does the wake begin to develop 
streamwise structures similar in intensity to the wakes initiated with pairs of oblique 
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Figure G.l: Iso-enstrophy density contour for 5° oblique fundamental plus two-dim- 
ensional subharmonic. Rc{, = 180. |u?| = 0.2. t = 106.6. 
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Figure G.2: Iso-enstrophy density contour for 15° oblique fundamental plus two-dim- 
ensional subharmonic. /?e<, = 180. |w| = 0.2. t = 95.2. 
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Figure G.3: Iso-enstrophy density contour for 30° oblique fundamental plus two-dim- 
ensional subharmonic. Rtb = 180. Ul = 0.2. t = 110.5. 
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Figure G.4: Spanwise vorticity at various spanwise locations for 5° oblique funda- 
mental plus two-dimensional subharmonic. Re b = 180. t = 106.6. Contours are 
0.01 < |u> 3 | < 0.4 in increments of 0.05. 



30° oblique fundamental 
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Figure G.7: Iso-enstrophy density contour for 30° oblique fundamental plus two-dim- 
ensional subharmonic. Rtt, — 1800. |u.'| = 0.8. t = 101.9. 
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subharmonic waves at the same development time (compare the wake in figures G.l, 
G.2, and G.3 to the the Re b = 119 wake in figure 4.13a in section 4.3.2 on page 109). 

This result is consistent with the stability analysis of Flemming [16], which predicts 
strong growth of subharmonic disturbances at angles between approximately 45° and 
70°. The 5° oblique fundamental case produces oblique subharmonic disturbances 
at angles well below 45°, and thus the three-dimensional structures do not grow 
in strength. At an oblique shedding angle of 30° however, the resulting oblique 
subharmonic disturbances are over 45° (see figure G.5) and thus grow in strength. 

This is very significant to the development of fine scale motions at higher Reynolds 
numbers. Figures G.6 and G.l show iso-enstrophy density contours for the 5 and 
30° oblique fundamental wakes, respectively, but this time at a Reynolds number of 
Rc b = 1800. This gives a Reynolds number and development time roughly equivalent 
to the wake in figure 4.13d in section 4.3.2 on page 110. At the time shown, the 
5° wake, which does not produce significant streamwise structures, has developed 
no detectable fine scale motions even at this rather high Reynolds number. The 
only effect of increasing the Reynolds number has been to slow the diffusion of the 
vorticity in the spanwise rollers. In contrast, the 30° wake, which does produce strong 
coherent streamwise structures, develops strong fine scale motions. This is in line with 
the earlier simulations, and reinforces the result that the tertiary transition in the 
incompressible plane wake requires the presence of strong long wavelength coherent 
streamwise structures. These results indicate that the source and symmetry of those 
structures are not significant, so long as they are present. 
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£ooi 

<^001 

£ 100 
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- 
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- 

- 

- 
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35 

283.4 

8 x 128 x 4 

- 
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- 

- 

- 
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- 
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69 
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Table H.l: Summary of Two-dimensional Simulations 

(Cont.) 


Tag 

Filename 


Ret, 

Max Time 
Ni xN 2 xN 3 


£100 


£010 

<^010 


£001 
4 > 001 


-100 


-010 


-001 


4> 


100 


4>' 


,010 


<f> { 


,001 


xxx 

Qxx 


119(0) 
wk2d_3 


119 

361.0 

128 x 128 x 4 


0.02 


xxx 

Qxx 


692(0) 
wk2d_4 


692 

337.8 

128 x 192 x 4 


0.02 


1384(0); 

wk2d_5 


XXX 

Qxx 


1384 

303.9 

192 x 256 x 4 


0.02 


2768(0); 

wk2d_6 


XXX 

Oxx 


2768 

291.0 

192 x 256 x 4 


0.02 


XXX 

OOx 


346(0) 
wk2d_7 


346 ( 0 ) 5 |* 

wk2d_8 


346 

311.4 

192 x 192 x 4 


0.02 


0.02 

0 


346 

251.5 

192 x 192 x 4 


0.02 


0.02 
7T / 4 


346 ( 0 ) 5 |^ 


)nr 

wk2d-9 


346 

312.4 

192 x 192 x 4 


0.02 


0.02 

tt/2 


XXX 

OOx 


119(0) 
wk2d-10 


119 

297.9 

128 x 128 x 4 


0.02 


0.02 

0 


















APPENDIX H. SUMMARY OF SIMULATIONS 


280 


Table H.l: Summary of Two-dimensional Simulations 

(Cont.) 
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Filename 
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Table H.l: Summary of Two-dimensional Simulations 
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Table H.2: Summary of Three-dimensional Simulations 
(Cont.) 
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Table H.2: Summary of Three-dimensional Simulations 
(Cont.) 
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Table H.2: Summary of Three-dimensional Simulations 
(Cont.) 
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Table H.2: Summary of Three-dimensional Simulations 
(Cont.) 
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Table H.2: Summary of Three-dimensional Simulations 
(Cont.) 
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